********************************************************************************
* 
********************************************************************************
clear programs
clear all
set more off

capture program drop main
program define main
	disp "Calling subprograms"
	PrepRDData
	RD_Plots_MainPaper
	RD_Plots_BilledElectricity
	RD_Plots_BalanceTests
	RD_FirstStage2ndStage
	RD_1stStage2ndStage_Chal23
	RD_McCrary
	RD_Ests
	iv_lnkwhw

end

*===============================================================================
capture program drop PrepRDData
program define PrepRDData
*===============================================================================

cd "$GeneralModified"
use "BalancedData",replace

*-------------------------------------------------------------------------------
* Simplify dataset for RD estimates
*-------------------------------------------------------------------------------
*drop indicators and merge in the total_kwh_num and similar values
capture drop t_t*
foreach x of numlist 0(1)10 {
	capture drop t_b`x'_fm*
	capture drop t_a`x'_fm*
}
foreach x of numlist 4(1)10 {
	capture drop t_b`x'*
	capture drop t_a`x'*
}

drop reduc_gap_*
drop lastmonth
drop reward

*Merge in info needed like total_kwh_num
rename reduc challenge_order
merge m:1 id challenge_order using Accounts_All,keepusing(total_kwh_num adjusted_historical_kwh_num adjusted_historical_kwh_num success2) keep(1 3) nogen
rename challenge_order reduc

*-------------------------------------------------------------------------------
* Indicator for whether HH continues past the current id or not
*-------------------------------------------------------------------------------
gen Continue=1 if num_challenges>reduc & reduc!=0
replace Continue=0 if Continue==.

bysort id reduc: gen uniid_reduc=1 if _n==1

*In each challenge and post challenge year 1-4, what is the total use as judged from monthly?
gen prepostyears=0 if t_a0==1
replace prepostyears=-3 if t_b3==1
replace prepostyears=-2 if t_b2==1
replace prepostyears=-1 if t_b1==1
replace prepostyears=1 if t_a1==1
replace prepostyears=2 if t_a2==1
replace prepostyears=3 if t_a3==1

bysort id prepostyears: gen uniprepost=1 if _n==1

gen numobs=.
bysort id prepostyears: replace numobs=_N if prepostyears!=.

*===============================================================================
* Compare BCHsuccess and BCHchange
*===============================================================================
*Generate what BCH evaluates the % change is. This should match success2 as judged via the 0.095% threshold
gen BCHchange=(total_kwh_num-adjusted_historical_kwh_num)/adjusted_historical_kwh_num if reduc!=0
gen BCHevalSuccess=1 if BCHchange<-0.095
replace BCHevalSuccess=0 if BCHevalSuccess==.
// tab BCHevalSuccess success2

*===============================================================================
* Save the BCH_2_evalSuccess measure of success in the 2nd challenge to go with 1st challenge
*===============================================================================
tab success2 BCHevalSuccess if reduc==1 & uniid_reduc==1
*Generate a new BCH_2_evalSuccess measure that applies across all observations not just during the given reduction year
capture drop temp
gen temp=1 if success2==1 & reduc==2
replace temp=0 if temp==. & reduc==2
capture drop BCH_2_evalSuccess
bysort id: egen BCH_2_evalSuccess=mean(temp)
// su BCH_2_evalSuccess
// tab success2 BCH_2_evalSuccess if reduc==2 & uniid_reduc==1

*-------------------------------------------------------------------------------
*Gen Centered BChange
*-------------------------------------------------------------------------------
gen Center_BCHchange=BCHchange+0.095
//what was -0.095 is now defined as 0

gen Center_BCHchange_s1=Center_BCHchange if BCHevalSuccess==1
replace Center_BCHchange_s1=0 if Center_BCHchange_s1==.
*tw sc Center_BCHchange_s1 Center_BCHchange if reduc==1 & uniid_reduc==1
*tw sc Center_BCHchange_s1 Center_BCHchange if reduc==1 & uniid_reduc==1 & Center_BCHchange<2,xline(0)

*-------------------------------------------------------------------------------
*Extend the Center_BCHchange to other observations beyond reduc==1
*-------------------------------------------------------------------------------
*Center_BCHchange is for the 1st challenge
capture drop temp
capture drop id_Center_BCHchange
bysort id: egen temp=mean(Center_BCHchange) if reduc==1
bysort id: egen id_Center_BCHchange=mean(temp)
capture drop temp

*===============================================================================
* Generate mychange and y2y1y0pcchange
*===============================================================================
* y2y1y0pcchange the change in use to year a1 (2nd challenge) from the 12 months prior to this.
* y2y1y0pcchange IS NOT the post-first challenge change!
* y2y1y0pcchange_v2 is defined as the 12 months post first challenge minus the 1st challenge use, normalized by t_b1 

*Need to ensure these totals are only calculated when there are 12 months of observations
bysort id prepostyears: gen num_months=_N if prepostyears!=.

*bysort id prepostyears: gen uniprepost=1 if _n==1
*IS THIS SECTION CORRECT TO USE shortgap==1 instead of months==1?
//likely yes; here, I"m interested in the decision to continue (within a short time period) and comparing the t_a0 to next months use.
//need to verify results are continuing within 12 months. NO.
//here I'm just defining the change in kwhw from one challenge to the following year. Later, in the regressions, I use shortgapV2 to define the regression only for those re-enrolling within 12 months. 

capture drop temp
forvalues i=0/3 {
	disp "`i'"
	capture drop temp
	bysort id: egen temp=total(kwhw) if t_a`i'==1 & num_months==12
	bysort id: egen ta`i'_year_kwhw_1=mean(temp)
}
forvalues i=1/2 {
	disp "`i'"
	capture drop temp
	bysort id: egen temp=total(kwhw) if t_b`i'==1 & num_months==12
	bysort id: egen tb`i'_year_kwhw_1=mean(temp)
}

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
* New way to calculate use in the 12 months prior to t_a1 and define tot_bChal2
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

*restore,not
*preserve

capture drop firstmonth
bysort id reduc (date): gen firstmonth=1 if _n==1 & reduc==2

tsset id date
capture drop t_t_bfa1
gen t_t_bfa1=.
forvalues i=1/12 {
	replace t_t_bfa1=1 if F`i'.firstmonth==1
}

capture drop temp2
capture drop tot_bChal2
bysort id: egen temp2=total(kwhw) if t_t_bfa1==1
replace temp2=ta0_year_kwhw_1 if reduc==1 & num_challenges==1
bysort id: egen tot_bChal2=mean(temp2)
capture drop temp2

capture drop t_t_bfa1
capture drop firstmonth

*Calculate the 12 months after the initial challenge
capture drop t_t_afa0
gen t_t_afa0=0
replace t_t_afa0=1 if L12.t_a0==1

capture drop temp2
capture drop tot_aChal1
bysort id: egen temp2=total(kwhw) if t_t_afa0==1
bysort id: egen tot_aChal1=mean(temp2)
capture drop temp2
capture drop t_t_afa0
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
* Generate y2y1y0pcchange using the t_t_b totals
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*br id tps date reduc num_challenges t_b* prepost* BCHchange tb* ta* t_t* tot* if tps==1 & num_challenges>=2
*br id tps date kwhw reduc num_challenges t_b* t_a* prepost* BCHchange tb* ta* t_t* tot* y2* if tps==1 & num_challenges>=2
*id 20004

*Need to make this the prior 12 months substracted from challegne 2, not the 12 months of challenge 1
*That is, find the total for ta1 lagged 12 months. This differes based on whether somebody continues to a second challenge or not
gen y2y1y0pcchange=(ta1_year_kwhw_1-tot_bChal2)/(tb1_year_kwhw_1) if reduc==1
*gen y2y1y1pcchange=(ta1_year_kwhw_1-tot_bChal2)/(tot_bChal2) if reduc==1
//why not define the change w.r.t.y1 as the baseline, not y0? 

*calculate the chagne from t_a0 to the 12 months after this, regardless of whether the household undertakes another conservation challenge
gen y2y1y0pcchange_v2=(tot_aChal1-ta0_year_kwhw_1)/(tb1_year_kwhw_1) if reduc==1
*tw sc y2y1y0pcchange_v2 y2y1y0pcchange if uniid_reduc==1 & reduc==1

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
* Generate y3y2y0pcchange
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*calculate the change TO t_a2 (to the third challenge) from the 12 months prior to this
tsset
capture drop t_t_bfa2
gen t_t_bfa2=.
replace t_t_bfa2=1 if F12.t_a2==1

capture drop temp2
capture drop tot_bChal3
bysort id: egen temp2=total(kwhw) if t_t_bfa2==1
bysort id: egen tot_bChal3=mean(temp2)
capture drop temp2

*check if have 12 months of t_a2 (third challenge) in which to calculate use
bysort id: egen tot_t_a2=total(t_a2)

gen y3y2y0pcchange=(ta2_year_kwhw_1-tot_bChal3)/(tb1_year_kwhw_1) if reduc==2 & tot_t_a2==12
capture drop t_t_bfa2

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
* Generate y4y3y0pcchange
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
tsset
capture drop t_t_bfa3
gen t_t_bfa3=.
replace t_t_bfa3=1 if F12.t_a3==1

capture drop temp2
capture drop tot_bChal4
bysort id: egen temp2=total(kwhw) if t_t_bfa3==1
bysort id: egen tot_bChal4=mean(temp2)
capture drop temp2

bysort id: egen tot_t_a3=total(t_a3)

gen y4y3y0pcchange=(ta3_year_kwhw_1-tot_bChal4)/(tb1_year_kwhw_1) if reduc==3 & tot_t_a3==12
capture drop t_t_bfa3

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
* Generate mychange
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
gen mychange=(ta0_year_kwhw_1-tb1_year_kwhw_1)/(tb1_year_kwhw_1) if reduc==1
bysort id: egen id_mychange=mean(mychange)

gen Center_mychange=mychange+0.095


*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
* Define the pre-program use and start_hdd to check continuity in covariates at the discontinuity
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
capture drop temp
bysort id t_b1: egen temp=total(kwhw) if t_b1==1
bysort id: egen t_b1totaluse=mean(temp)
drop temp

capture drop temp
bysort id t_b2: egen temp=total(kwhw) if t_b2==1
bysort id: egen t_b2totaluse=mean(temp)
drop temp

cd "$datafrom"
rename date year_month
merge m:1 year_month using hdd_data,nogen keep(3)
rename year_month date

capture drop temp
bysort id: gen temp=hdd if date==begin_challenge_1-1
capture drop start_hdd
bysort id: egen start_hdd=mean(temp)

capture drop temp
bysort id: egen temp=total(hdd) if t_b1==1
bysort id: egen start_yr_hdd=mean(temp)

capture drop temp
bysort id: egen temp=total(hdd) if t_a0==1
bysort id: egen start_yr_hdda0=mean(temp)

drop hdd
drop temp

*For month
gen start_month=month(dofm(begin_challenge_1))

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
* Simplify to cross-section data
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
drop tot_t_a2
drop tot_t_a3
keep if tps==1 & uniid_reduc==1

cd "$GeneralModified"
save "temp_RD_crosssection",replace

end

*===============================================================================
capture program drop RD_Plots_MainPaper
program define RD_Plots_MainPaper
*===============================================================================
*Bin width 0.0075. Includes both FS and RF plots

cd "$GeneralModified"
use "temp_RD_crosssection",replace
*-------------------------------------------------------------------------------
* Generate fixed width bins of BCHchange — Defines the x axis
*-------------------------------------------------------------------------------
*Define a bin width and let the # bins correspond to the binwidth. 
capture drop bins_BCHfixedbins
gen bins_BCHfixedbins=.
local binwidth=0.0075
local i=1
forvalues x=-0.095(-`binwidth')-0.30 {
	local i=`i'-1
	local x1=`x'-`binwidth'
	replace bins_BCHfixedbins=`i' if inrange(BCHchange,`x1',`x')
}
*Find the smallest BCHfixedbin value and use this to offset bin numbers when counting above -0.095
su bins_BCHfixedbins
local minval=`r(min)'
disp "`minval'"
replace bins_BCHfixedbins=bins_BCHfixedbins-`minval'+1 if bins_BCHfixedbins!=.
su bins_BCHfixedbins
local i=`r(max)'
forvalues x=-0.095(`binwidth')0.2 {
	local i=`i'+1
	local x1=`x'+`binwidth'
	replace bins_BCHfixedbins=`i' if inrange(BCHchange,`x',`x1')
}

*-------------------------------------------------------------------------------
*Save the mean values of BCHchange by bin number, bins_BCHfixedbins.
//bins_BCHfixedbins are ordered bin numbers X_BCHfixedbins1
*-------------------------------------------------------------------------------
su bins_BCHfixedbins
local maxval=`r(max)'
sort bins_BCHfixedbins

forvalues chal=1/1 {
	capture matrix drop A
	capture matrix drop temp
	tabstat BCHchange if uniid_reduc==1 & reduc==`chal' & balancedset_sg6==1,by(bins_BCHfixedbins) save
	matrix A=.
	foreach mval of numlist 1/`maxval' {
		matrix temp=r(Stat`mval')
		matrix A=A\temp
	}
	matrix A=A[2...,1]
	capture drop X_BCHfixedbins`chal'
	svmat A,names(X_BCHfixedbins`chal')
	matrix X_BCHfixedbins`chal'=A
	rename X_BCHfixedbins`chal'1 X_BCHfixedbins`chal'
}

capture drop ests_discont_BCH_*
capture drop CIupper_dc_BCH_* CIlower_dc_BCH_*
capture drop stder_dc_BCH_*

*-------------------------------------------------------------------------------
*Plots: average of Continue by bin number
*-------------------------------------------------------------------------------
sort X_BCHfixedbins1
eststo discont_1: quietly reg Continue ibn.bins_BCHfixedbins if uniid_reduc==1 & reduc==1 & balancedset==1 & shortgapV2<=12,nocon vce(robust)
sort X_BCHfixedbins1
matrix ests=e(b)'
svmat ests,names(ests_discont_BCH_1)
rename ests_discont_BCH_11 ests_discont_BCH_1

matrix ses=e(V)'
matrix ses=vecdiag(ses)'
svmat ses,names(stder_dc_BCH_1)
rename stder_dc_BCH_11 stder_dc_BCH_1
replace stder_dc_BCH_1=sqrt(stder_dc_BCH_1)

gen CIupper_dc_BCH_1=ests_discont_BCH_1+stder_dc_BCH_1*1.96
gen CIlower_dc_BCH_1=ests_discont_BCH_1-stder_dc_BCH_1*1.96

capture drop rarea*

gen rarea_x=-0.093 if _n==1
replace rarea_x=0.2 if _n==2
gen rarea_y1=1.03 if _n<=2
gen rarea_y2=1.1 if _n<=2

gen rarea_x2=-0.097 if _n==1
replace rarea_x2=-0.3 if _n==2
gen rarea_y12=1.03 if _n<=2
gen rarea_y22=1.1 if _n<=2

tw rspike CIupper_dc_BCH_1 CIlower_dc_BCH_1 X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<0.205,lstyle(ci) lcolor(navy) yscale(range(0 1.08) axis(1)) ///
	|| sc ests_discont_BCH_1 X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<0.205, mcolor(navy)  ///
	|| lpoly ests_discont_BCH_1 X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<-0.095,deg(1) lp(dash) lc(gs2) lw(medthin) mc(gs4) msize(vsmall) yscale(range(0 1)) ///
	|| lpoly ests_discont_BCH_1 X_BCHfixedbins1 if X_BCHfixedbins1>-0.095 & X_BCHfixedbins1<0.205,deg(1) lp(dash) lc(gs2) lw(medthin) mc(gs4) msize(vsmall) yscale(range(0 1)) ///
	|| rarea rarea_y1 rarea_y2 rarea_x, sort color(red*0.2) ///
	|| rarea rarea_y12 rarea_y22 rarea_x2, sort color(green*0.2) ///
	ylabel(0(0.2)1,axis(1) glcolor(gs10)) ///
	ytitle("Probability of Re-enrolling",axis(1)) ///
	xlabel(none) xtitle(Credited Changes) ///
	xlabel(-0.3(0.1).2 -0.3 "-30%" -0.2 "-20%" -0.1 "-10%" 0 "0%" 0.1 "10%" 0.2 "20%") ///
	xtick(-0.3(0.1).2) xline(-0.095,lc(maroon) lp(dash)) graphregion(color(white)) legend(off) ///
	text(1.07 -0.2 "Pass 1st Challenge") text(1.07 0.05 "Fail 1st Challenge") 
	
graph rename Simple_S1_b075_area_v4,replace
//NOT IN PAPER — full version of first stage plot

*-------------------------------------------------------------------------------
*Combine the first stage and RF plots into a single figure, save this FStemp plot for later combining
gen rarea_y1_win=1.01
gen rarea_y2_win=1.08
gen rarea_y12_win=1.01
gen rarea_y22_win=1.08
	
tw rspike CIupper_dc_BCH_1 CIlower_dc_BCH_1 X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<0.205,lstyle(ci) lcolor(navy) yscale(range(0 1) axis(1)) ///
	|| sc ests_discont_BCH_1 X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<0.205, mcolor(navy)  ///
	|| lpoly ests_discont_BCH_1 X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<-0.095,deg(1) lp(dash) lc(gs2) lw(medthin) mc(gs4) msize(vsmall) yscale(range(0 1)) ///
	|| lpoly ests_discont_BCH_1 X_BCHfixedbins1 if X_BCHfixedbins1>-0.095 & X_BCHfixedbins1<0.205,deg(1) lp(dash) lc(gs2) lw(medthin) mc(gs4) msize(vsmall) yscale(range(0 1)) ///
	|| rarea rarea_y1_win rarea_y2_win rarea_x, sort color(red*0.2) ///
	|| rarea rarea_y12_win rarea_y22_win rarea_x2, sort color(green*0.2) ///
	ylabel(0(0.2)1,axis(1) glcolor(gs10)) ///
	ytitle("Probability of Re-enrolling",axis(1)) ///
	text(1.05 -0.2 "Pass 1st Challenge") text(1.05 0.05 "Fail 1st Challenge")  fysize(60) text(1.13 -.07  "(a) First Stage",size(medsmall)) ///
	xtick(-0.3(0.1).2) xline(-0.095,lc(maroon) lp(dash)) graphregion(color(white)) legend(off) ///
	xlabel(none) xtitle("") plotr(margin(l r b t-2))
graph rename FStemp,replace		
//COMBINED BELOW AND IS IN PAPER

*-------------------------------------------------------------------------------
*Plots - Second challenge discontinuity in continuing
*-------------------------------------------------------------------------------

sort X_BCHfixedbins1
eststo discont_2: quietly reg Continue ibn.bins_BCHfixedbins if uniid_reduc==1 & reduc==2 & balancedset==1 & shortgapV2<=12,nocon vce(robust)
sort X_BCHfixedbins1
matrix ests=e(b)'
capture drop ests_discont_BCH_2*
capture drop _ests_discont_BCH_2*
svmat ests,names(ests_discont_BCH_2)
rename ests_discont_BCH_21 ests_discont_BCH_2

matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stder_dc_BCH_2*
svmat ses,names(stder_dc_BCH_2)

rename stder_dc_BCH_21 stder_dc_BCH_2
replace stder_dc_BCH_2=sqrt(stder_dc_BCH_2)

capture drop CIupper_dc_BCH_2*
capture drop CIlower_dc_BCH_2*
gen CIupper_dc_BCH_2=ests_discont_BCH_2+stder_dc_BCH_2*1.96
gen CIlower_dc_BCH_2=ests_discont_BCH_2-stder_dc_BCH_2*1.96

tw rspike CIupper_dc_BCH_2 CIlower_dc_BCH_2 X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<0.205,lstyle(ci) lcolor(navy) yscale(range(0 1.08) axis(1)) ///
	|| sc ests_discont_BCH_2 X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<0.205, mcolor(navy)  ///
	|| lpoly ests_discont_BCH_2 X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<-0.095,deg(1) lp(dash) lc(gs2) lw(medthin) mc(gs4) msize(vsmall) yscale(range(0 1)) ///
	|| lpoly ests_discont_BCH_2 X_BCHfixedbins1 if X_BCHfixedbins1>-0.095 & X_BCHfixedbins1<0.205,deg(1) lp(dash) lc(gs2) lw(medthin) mc(gs4) msize(vsmall) yscale(range(0 1)) ///
	ylabel(0(0.2)1,axis(1) glcolor(gs10)) ///
	ytitle("Probability of Re-enrolling",axis(1)) ///
	xlabel(none) xtitle(Credited Changes) ///
	xlabel(-0.3(0.1).2 -0.3 "-30%" -0.2 "-20%" -0.1 "-10%" 0 "0%" 0.1 "10%" 0.2 "20%") ///
	xtick(-0.3(0.1).2) xline(-0.095,lc(maroon) lp(dash)) graphregion(color(white)) legend(off) ///
	text(1.15 -0.2 "Pass 2nd Challenge") text(1.15 0.05 "Fail 2nd Challenge") 

graph rename RD_Continue_2ndChallenge,replace
cd "$SaveOutputFiles"
graph export RD_Continue_2ndChallenge.png,replace
//Figure F.1 (a)

*-------------------------------------------------------------------------------
*Plots - Third challenge discontinuity in continuing
*-------------------------------------------------------------------------------

sort X_BCHfixedbins1
eststo discont_3: quietly reg Continue ibn.bins_BCHfixedbins if uniid_reduc==1 & reduc==3 & balancedset==1 & shortgapV2<=12,nocon vce(robust)
sort X_BCHfixedbins1
matrix ests=e(b)'
capture drop ests_discont_BCH_3*
capture drop _ests_discont_BCH_3*
svmat ests,names(ests_discont_BCH_3)
rename ests_discont_BCH_31 ests_discont_BCH_3

matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stder_dc_BCH_3*
svmat ses,names(stder_dc_BCH_3)
rename stder_dc_BCH_31 stder_dc_BCH_3
replace stder_dc_BCH_3=sqrt(stder_dc_BCH_3)

capture drop CIupper_dc_BCH_3*
capture drop CIlower_dc_BCH_3*
gen CIupper_dc_BCH_3=ests_discont_BCH_3+stder_dc_BCH_3*1.96
gen CIlower_dc_BCH_3=ests_discont_BCH_3-stder_dc_BCH_3*1.96

tw rspike CIupper_dc_BCH_3 CIlower_dc_BCH_3 X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<0.205,lstyle(ci) lcolor(navy) yscale(range(0 1.08) axis(1)) ///
	|| sc ests_discont_BCH_3 X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<0.205, mcolor(navy)  ///
	|| lpoly ests_discont_BCH_3 X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<-0.095,deg(1) lp(dash) lc(gs2) lw(medthin) mc(gs4) msize(vsmall) yscale(range(0 1)) ///
	|| lpoly ests_discont_BCH_3 X_BCHfixedbins1 if X_BCHfixedbins1>-0.095 & X_BCHfixedbins1<0.205,deg(1) lp(dash) lc(gs2) lw(medthin) mc(gs4) msize(vsmall) yscale(range(0 1)) ///
	ylabel(0(0.2)1,axis(1) glcolor(gs10)) ///
	ytitle("Probability of Re-enrolling",axis(1)) ///
	xlabel(none) xtitle(Credited Changes) ///
	xlabel(-0.3(0.1).2 -0.3 "-30%" -0.2 "-20%" -0.1 "-10%" 0 "0%" 0.1 "10%" 0.2 "20%") ///
	xtick(-0.3(0.1).2) xline(-0.095,lc(maroon) lp(dash)) graphregion(color(white)) legend(off) ///
	text(1.2 -0.2 "Pass 3rd Challenge") text(1.2 0.05 "Fail 3rd Challenge") 
	
graph rename RD_Continue_3rdChallenge,replace
cd "$SaveOutputFiles"
graph export RD_Continue_3rdChallenge.png,replace
//Figure F.1 (b)

*-------------------------------------------------------------------------------
*Estimates - Post-challenge change 
*-------------------------------------------------------------------------------
eststo m13_1_both: quietly reg y2y1y0pcchange ibn.bins_BCHfixedbins if uniid_reduc==1 & reduc==1 & balancedset==1 & shortgapV2<=12,nocon vce(robust)

sort X_BCHfixedbins1
capture drop ests_m13_1_both*
matrix ests=e(b)'
svmat ests,names(ests_m13_1_both)

matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_m13_1_both*
capture drop CIupper_m13_1_both*
capture drop CIlower_m13_1_both*
svmat ses,names(stderrors_m13_1_both)
replace stderrors_m13_1_both=sqrt(stderrors_m13_1_both)
gen CIupper_m13_1_both=ests_m13_1_both+stderrors_m13_1_both*1.96
gen CIlower_m13_1_both=ests_m13_1_both-stderrors_m13_1_both*1.96

*May need to winsorize the data for plotting it to make it narrower
gen CIupper_m13_1_both_win=CIupper_m13_1_both
replace CIupper_m13_1_both_win=0.08 if CIupper_m13_1_both>=0.08
 
gen CIlower_m13_1_both_win=CIlower_m13_1_both
replace CIlower_m13_1_both_win=-0.10 if CIlower_m13_1_both<=-0.10

tw sc ests_m13_1_both X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<0.205 & ests_m13_1_both>-0.1,mcolor(navy)  /// 
 || rspike CIupper_m13_1_both_win CIlower_m13_1_both_win X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<0.205,lstyle(ci) lcolor(navy) ///
	|| lpoly ests_m13_1_both X_BCHfixedbins1 if X_BCHfixedbins1>-0.305 & X_BCHfixedbins1<-0.095,deg(1) lp(dash) lc(gs2) lw(medthin) mc(gs4) msize(vsmall) ///
	|| lpoly ests_m13_1_both X_BCHfixedbins1 if X_BCHfixedbins1>-0.095 & X_BCHfixedbins1<0.205,deg(1) lp(dash) lc(gs2) lw(medthin) mc(gs4) msize(vsmall) ///
	xlabel(-0.3(0.1).2 -0.3 "-30%" -0.2 "-20%" -0.1 "-10%" 0 "0%" 0.1 "10%" 0.2 "20%") ///
	xline(-0.095,lc(maroon) lp(dash)) graphregion(color(white)) legend(off) ///
	ylabel(-0.075(0.05)0.08  -0.075 "-7.5%" -0.025 "-2.5%" 0.075 "7.5%" 0.025 "2.5%" ,axis(1)) ///
	xtitle("Credited Changes during the inital challenge") ytitle("Post-challenge percent change in kWh") legend(off) graphregion(color(white)) ///
	yscale(range(-0.1 0.075)) fysize(40) text(.1 -.07  "(b) Reduced Form",size(medsmall))
graph rename RFtemp,replace

graph combine FStemp RFtemp,r(2) graphregion(color(white)) iscale(0.6) imargin(vsmall) name(RD_FS_RF,replace)
graph display RD_FS_RF, xsize(9) ysize(10)
//Paper Figure 2

end

*===============================================================================
capture program drop RD_Plots_BilledElectricity
program define RD_Plots_BilledElectricity
*===============================================================================
*Uses billed, not credited, changes in use. So used bins_myfixedbins to define bins of annual changes in kWh as measured in the data.
cd "$GeneralModified"
use "temp_RD_crosssection",replace

*===============================================================================
* Generate fixed width bins of mychange
*===============================================================================

capture drop bins_myfixedbins
gen bins_myfixedbins=.
local binwidth=0.0075
local i=1
forvalues x=-0.095(-`binwidth')-0.30 {
	local i=`i'-1
	disp("X is `x'")
	disp("I is `i'")
	local x1=`x'-`binwidth'
	disp("X1 is `x1'")
	replace bins_myfixedbins=`i' if inrange(mychange,`x1',`x')
}
*Find the minimul value and use this to offset the set of bins
su bins_myfixedbins
local minval=`r(min)'
disp "`minval'"
replace bins_myfixedbins=bins_myfixedbins-`minval'+1 if bins_myfixedbins!=.
su bins_myfixedbins
local i=`r(max)'
forvalues x=-0.095(`binwidth')0.2 {
	local i=`i'+1
	disp("X is `x'")
	disp("I is `i'")
	local x1=`x'+`binwidth'
	replace bins_myfixedbins=`i' if inrange(mychange,`x',`x1')
}

su bins_myfixedbins if tps==1 & reduc==1

*-------------------------------------------------------------------------------
*Save the mean values of bins_myfixedbins
*-------------------------------------------------------------------------------
su bins_myfixedbins
local maxval=`r(max)'
sort bins_myfixedbins

forvalues chal=1/1 {
	capture matrix drop A
	capture matrix drop temp
	tabstat mychange if uniid_reduc==1 & reduc==`chal' & balancedset==1 & shortgapV2<=12,by(bins_myfixedbins) save
	matrix A=.
	foreach mval of numlist 1/`maxval' {
		matrix temp=r(Stat`mval')
		matrix A=A\temp
	}
	matrix A=A[2...,1]
	capture drop X_myfixedbins`chal'
	svmat A,names(X_myfixedbins`chal')
	matrix X_myfixedbins`chal'=A
	rename X_myfixedbins`chal'1 X_myfixedbins`chal'
}
	
*-------------------------------------------------------------------------------
* Plot vs mychange bins
*-------------------------------------------------------------------------------

capture drop ests_discont_my_*
capture drop CIupper_dc_my_* CIlower_dc_my_*
capture drop stder_dc_my_*

*Record the P(continue) from a regression for equal width bins
sort X_myfixedbins1
eststo discont_1: quietly reg Continue ibn.bins_myfixedbins if uniid_reduc==1 & reduc==1 & balancedset==1 & shortgapV2<=12,nocon vce(robust)

matrix ests=e(b)'
svmat ests,names(ests_discont_my_1)
rename ests_discont_my_11 ests_discont_my_1

matrix ses=e(V)'
matrix ses=vecdiag(ses)'
svmat ses,names(stder_dc_my_1)
rename stder_dc_my_11 stder_dc_my_1
replace stder_dc_my_1=sqrt(stder_dc_my_1)

gen CIupper_dc_my_1=ests_discont_my_1+stder_dc_my_1*1.96
gen CIlower_dc_my_1=ests_discont_my_1-stder_dc_my_1*1.96

*For success
eststo m14_4_S1: quietly reg BCHevalSuccess ibn.bins_myfixedbins if uniid_reduc==1 & reduc==1 & balancedset==1 & shortgapV2<=12,nocon vce(robust)

sort X_myfixedbins1
capture drop ests_m14_4_S1*
matrix ests=e(b)'
svmat ests,names(ests_m14_4_S1)

matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_m14_4_S1*
capture drop CIupper_m14_4_S1*
capture drop CIlower_m14_4_S1*
svmat ses,names(stderrors_m14_4_S1)
replace stderrors_m14_4_S1=sqrt(stderrors_m14_4_S1)
gen CIupper_m14_4_S1=ests_m14_4_S1+stderrors_m14_4_S1*1.96
gen CIlower_m14_4_S1=ests_m14_4_S1-stderrors_m14_4_S1*1.96

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
* Include a third line that is the P(success) in the following challenge?
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 	
*For success in next challenge? 
eststo m14_4_S2: quietly reg BCH_2_evalSuccess ibn.bins_myfixedbins if uniid_reduc==1 & reduc==1 & balancedset==1 & shortgapV2<=12,nocon vce(robust)

sort X_myfixedbins1
capture drop ests_m14_4_S2*
matrix ests=e(b)'
svmat ests,names(ests_m14_4_S2)

matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_m14_4_S2*
capture drop CIupper_m14_4_S2*
capture drop CIlower_m14_4_S2*
svmat ses,names(stderrors_m14_4_S2)
replace stderrors_m14_4_S2=sqrt(stderrors_m14_4_S2)
gen CIupper_m14_4_S2=ests_m14_4_S2+stderrors_m14_4_S2*1.96
gen CIlower_m14_4_S2=ests_m14_4_S2-stderrors_m14_4_S2*1.96

*Alternate ordering of the panels
tw sc ests_m14_4_S1 X_myfixedbins1 if X_myfixedbins1>-0.305 & X_myfixedbins1<0.205, mc(gs2) msize(small) yscale(range(0 1)) ///
	|| line ests_m14_4_S1 X_myfixedbins1 if X_myfixedbins1>-0.305 & X_myfixedbins1<0.205, lc(gs2) msize(small) yscale(range(0 1)) ///
	|| sc ests_m14_4_S2 X_myfixedbins1 if X_myfixedbins1>-0.305 & X_myfixedbins1<0.205, mc(gs8) msize(small) yscale(range(0 1)) ///
	|| lpoly ests_m14_4_S2 X_myfixedbins1 if X_myfixedbins1>-0.305 & X_myfixedbins1<0.205,deg(3) lp(dash) lc(gs2) lw(medthin) mc(gs4) msize(vsmall) yscale(range(0 1)) ///
	ylabel(0(0.25)1,axis(1)) xtitle("") ///
	ytitle(Fraction that Pass,axis(1)) xtitle(Billed Changes) ///
	xlabel(-0.3(0.1).2 -0.3 "-30%" -0.2 "-20%" -0.1 "-10%" 0 "0%" 0.1 "10%" 0.2 "20%") fysize(32) ///
	xline(-0.095,lc(maroon) lp(dash)) graphregion(color(white)) legend(order(1 "1st Challenge" 3 "2nd Challenge"))
graph rename Simple_S1_3_S2,replace

tw rspike CIupper_dc_my_1 CIlower_dc_my_1 X_myfixedbins1 if X_myfixedbins1>-0.305 & X_myfixedbins1<0.205,lstyle(ci) lcolor(navy) yscale(range(0 1.08) axis(1)) ///
	|| sc ests_discont_my_1 X_myfixedbins1 if X_myfixedbins1>-0.305 & X_myfixedbins1<0.205, mcolor(navy) ///
	|| lpoly ests_discont_my_1 X_myfixedbins1 if X_myfixedbins1>-0.305 & X_myfixedbins1<-0.095,deg(1) lp(dash) lc(gs2) lw(medthin) mc(gs4) msize(vsmall) yscale(range(0 1)) ///
	|| lpoly ests_discont_my_1 X_myfixedbins1 if X_myfixedbins1>-0.095 & X_myfixedbins1<0.205,deg(1) lp(dash) lc(gs2) lw(medthin) mc(gs4) msize(vsmall) yscale(range(0 1)) ///
	ylabel(0(0.2)1,axis(1))  xtitle("") xlabel(none) ///
	ytitle(Probability of 2nd Challenge,axis(1)) ///
	xline(-0.095,lc(maroon) lp(dash)) graphregion(color(white)) legend(off)
graph rename Simple_S1_3_S1,replace

grc1leg Simple_S1_3_S1 Simple_S1_3_S2,r(2) graphregion(color(white)) iscale(1) imargin(small) legendfrom(Simple_S1_3_S2)
graph rename Simple_S1_3_combined_bin075_v2,replace

cd "$SaveOutputFiles"
graph export Simple_S1_3_combined_bin075_v2.png,replace
//Appendix Figure D.1

end

*===============================================================================
capture program drop RD_Plots_BalanceTests
program define RD_Plots_BalanceTests
*===============================================================================
*Use binwidth of 0.0025
cd "$GeneralModified"
use "temp_RD_crosssection",replace

*Merge in building characteristics for balance tests
cd "$datafrom"
merge m:1 id using unique_buildingchars_class,keepusing(floorarea  value_tot) keep(1 3) nogen
	
*-------------------------------------------------------------------------------
* Change the binsize of bins_BCHfixedbins_v2 if necessary
*-------------------------------------------------------------------------------
capture drop bins_BCHfixedbins_v2
gen bins_BCHfixedbins_v2=.
local binwidth=0.0025
local i=1
forvalues x=-0.095(-`binwidth')-0.195 {
	local i=`i'-1
	disp("X is `x'")
	disp("I is `i'")
	local x1=`x'-`binwidth'
	disp("X1 is `x1'")
	replace bins_BCHfixedbins_v2=`i' if inrange(BCHchange,`x1',`x')
}
*Find the minimul value and use this to offset the set of bins
su bins_BCHfixedbins_v2
local minval=`r(min)'
disp "`minval'"
replace bins_BCHfixedbins_v2=bins_BCHfixedbins_v2-`minval'+1 if bins_BCHfixedbins_v2!=.
su bins_BCHfixedbins_v2
local i=`r(max)'
forvalues x=-0.095(`binwidth')0.005 {
	local i=`i'+1
	disp("X is `x'")
	disp("I is `i'")
	local x1=`x'+`binwidth'
	replace bins_BCHfixedbins_v2=`i' if inrange(BCHchange,`x',`x1')
}

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
* Save mean values of bins_BCHfixedbins_v2 in narrow window of threshold (+/-5%) or plus/minus 10
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
su bins_BCHfixedbins_v2
local maxval=`r(max)'
sort bins_BCHfixedbins_v2

capture matrix drop A
capture matrix drop temp
tabstat BCHchange if uniid_reduc==1 & reduc==1 & balancedset_sg6==1,by(bins_BCHfixedbins_v2) save
matrix A=.
foreach mval of numlist 1/`maxval' {
	matrix temp=r(Stat`mval')
	matrix A=A\temp
}
matrix A=A[2...,1]
capture drop X_BCHfixedbins_v2
svmat A,names(X_BCHfixedbins_v2)
matrix X_BCHfixedbins_v2=A
rename X_BCHfixedbins_v21 X_BCHfixedbins_v2

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*Generate the negative of X_BCHfixedbins_v2 to flip the x axis in the plots
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
capture drop X_BCHfixedbins_v2_neg
gen X_BCHfixedbins_v2_neg=-X_BCHfixedbins_v2

*-------------------------------------------------------------------------------
*Check heating
*-------------------------------------------------------------------------------
tabstat heating if uniid_reduc==1 & reduc==1 & balancedset_sg6==1,by(bins_BCHfixedbins_v2)

*Check heating via regression
gen shareht=1 if heating==1
replace shareht=0 if heating!=1
eststo balance_heating: quietly reg shareht ibn.bins_BCHfixedbins_v2 if uniid_reduc==1 & reduc==1 & balancedset_sg6==1,nocon vce(robust)

sort X_BCHfixedbins_v2
matrix ests=e(b)'
capture drop ests_heat_*
svmat ests,names(ests_heat_)

matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_heat_*
svmat ses,names(stderrors_heat_)
replace stderrors_heat_1=sqrt(stderrors_heat_1)
capture drop CIupper_heat_1 CIlower_heat_1
gen CIupper_heat_1=ests_heat_1+stderrors_heat_1*1.96
gen CIlower_heat_1=ests_heat_1-stderrors_heat_1*1.96

tw rspike CIupper_heat_1 CIlower_heat_1 X_BCHfixedbins_v2 if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005,lstyle(ci) lcolor(navy) ///
	|| sc ests_heat_1 X_BCHfixedbins_v2 if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005, mcolor(navy) lc(navy) ///
	ytitle(Share With Electric Heating) xtitle("") xline(-0.095,lc(maroon) lp(dash)) xline(-0.1,lc(maroon) lp(solid)) ///
	xlabel(-0.2(0.05)0 0 "0%" -0.05 "-5%" -0.1 "-10%" -0.15 "-15%" -0.2 "-20%") graphregion(color(white)) ///
 	legend(off) ylabel(#5) title("Electric Heating")
graph rename BalanceTest_Heating,replace
//Figure E.3 (a)

*-------------------------------------------------------------------------------
*Check buildingtype
*-------------------------------------------------------------------------------
gen bt_SFD=1 if buildingtype==1
replace bt_SFD=0 if buildingtype!=1

eststo balance_buildingtype: quietly reg bt_SFD ibn.bins_BCHfixedbins_v2 if uniid_reduc==1 & reduc==1 & balancedset_sg6==1,nocon vce(robust)

sort X_BCHfixedbins_v2
matrix ests=e(b)'
capture drop ests_BT_*
svmat ests,names(ests_BT_)

matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_BT_*
svmat ses,names(stderrors_BT_)
replace stderrors_BT_1=sqrt(stderrors_BT_1)
capture drop CIupper_BT_1 CIlower_BT_1
gen CIupper_BT_1=ests_BT_1+stderrors_BT_1*1.96
gen CIlower_BT_1=ests_BT_1-stderrors_BT_1*1.96

tw rspike CIupper_BT_1 CIlower_BT_1 X_BCHfixedbins_v2_neg if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005,lstyle(ci) lcolor(navy) xline(0.095,lc(maroon) lp(dash)) ///
	|| con ests_BT_1 X_BCHfixedbins_v2_neg if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005, mcolor(navy) lc(navy) ///
	ytitle(Share SFD) xtitle("") ///
	xlabel(0(0.05)0.2 0 "0%" 0.05 "5%" 0.1 "10%" 0.15 "15%" 0.2 "20%") graphregion(color(white)) ///
	legend(off) ylabel(#5)
graph rename BalanceTest_buildingtype,replace
//not in paper

*-------------------------------------------------------------------------------	
*Floorarea
*-------------------------------------------------------------------------------
eststo balance_floorarea: quietly reg floorarea ibn.bins_BCHfixedbins_v2 if uniid_reduc==1 & reduc==1 & balancedset_sg6==1,nocon vce(robust)
sort X_BCHfixedbins_v2
matrix ests=e(b)'
capture drop ests_floorarea_*
svmat ests,names(ests_floorarea_)

matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_floorarea_*
svmat ses,names(stderrors_floorarea_)
replace stderrors_floorarea_1=sqrt(stderrors_floorarea_1)
capture drop CIupper_floorarea_1 CIlower_floorarea_1
gen CIupper_floorarea_1=ests_floorarea_1+stderrors_floorarea_1*1.96
gen CIlower_floorarea_1=ests_floorarea_1-stderrors_floorarea_1*1.96

tw rspike CIupper_floorarea_1 CIlower_floorarea_1 X_BCHfixedbins_v2 if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005,lstyle(ci) lcolor(navy) ///
	|| sc ests_floorarea_1 X_BCHfixedbins_v2 if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005, mcolor(navy) lc(navy) ///
	ytitle(Mean Floor Space) xtitle("") xline(-0.095,lc(maroon) lp(dash)) xline(-0.1,lc(maroon) lp(solid)) ///
	xlabel(-0.2(0.05)0 0 "0%" -0.05 "-5%" -0.1 "-10%" -0.15 "-15%" -0.2 "-20%") graphregion(color(white)) ///
 	legend(off) ylabel(#5) title("Floor Area")
graph rename Balance_Floorarea,replace
//Figure E.3 (b)

*-------------------------------------------------------------------------------		
*Total value
*-------------------------------------------------------------------------------
gen value_tot2=value_tot/1000
eststo balance_value_tot2: quietly reg value_tot2 ibn.bins_BCHfixedbins_v2 if uniid_reduc==1 & reduc==1 & balancedset_sg6==1,nocon vce(robust)
sort X_BCHfixedbins_v2
matrix ests=e(b)'
capture drop ests_value_*
svmat ests,names(ests_value_)
matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_value_*
svmat ses,names(stderrors_value_)
replace stderrors_value_1=sqrt(stderrors_value_1)
capture drop CIupper_value_1 CIlower_value_1
gen CIupper_value_1=ests_value_1+stderrors_value_1*1.96
gen CIlower_value_1=ests_value_1-stderrors_value_1*1.96

tw rspike CIupper_value_1 CIlower_value_1 X_BCHfixedbins_v2_neg if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005,lstyle(ci) lcolor(navy) xline(0.095,lc(maroon) lp(dash)) ///
	|| con ests_value_1 X_BCHfixedbins_v2_neg if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005, mcolor(navy) lc(navy) ///
	ytitle(Mean Household Value (Thousands)) xtitle("") ///
	xlabel(0(0.05)0.2 0 "0%" 0.05 "5%" 0.1 "10%" 0.15 "15%" 0.2 "20%") graphregion(color(white)) ///
	legend(off) ylabel(#5)
graph rename Balance_TotalValue,replace
//not in paper

*-------------------------------------------------------------------------------
* Use pre-program
*-------------------------------------------------------------------------------

eststo balance_tb2kwhw: quietly reg t_b2totaluse ibn.bins_BCHfixedbins_v2 if uniid_reduc==1 & reduc==1 & balancedset_sg6==1,nocon vce(robust)
sort X_BCHfixedbins_v2
matrix ests=e(b)'
capture drop ests_tb2kwhw_*
svmat ests,names(ests_tb2kwhw_)
matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_tb2kwhw_*
svmat ses,names(stderrors_tb2kwhw_)
replace stderrors_tb2kwhw_1=sqrt(stderrors_tb2kwhw_1)
capture drop CIupper_tb2kwhw_1 CIlower_tb2kwhw_1
gen CIupper_tb2kwhw_1=ests_tb2kwhw_1+stderrors_tb2kwhw_1*1.96
gen CIlower_tb2kwhw_1=ests_tb2kwhw_1-stderrors_tb2kwhw_1*1.96

tw rspike CIupper_tb2kwhw_1 CIlower_tb2kwhw_1 X_BCHfixedbins_v2_neg if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005,lstyle(ci) lcolor(navy) xline(0.095,lc(maroon) lp(dash)) ///
	|| con ests_tb2kwhw_1 X_BCHfixedbins_v2_neg if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005, mcolor(navy) lc(navy) ///
	ytitle(Mean kWh Pre-Challenge) xtitle("") ///
	xlabel(0(0.05)0.2 0 "0%" 0.05 "5%" 0.1 "10%" 0.15 "15%" 0.2 "20%") graphregion(color(white)) ///
	legend(off) ylabel(#5)
graph rename Balance_tb2kwhw,replace
//not in paper

*-------------------------------------------------------------------------------
* Use pre-program t_b1
*-------------------------------------------------------------------------------

eststo balance_tb1kwhw: quietly reg t_b1totaluse ibn.bins_BCHfixedbins_v2 if uniid_reduc==1 & reduc==1 & balancedset_sg6==1,nocon vce(robust)
sort X_BCHfixedbins_v2
matrix ests=e(b)'
capture drop ests_tb1kwhw_*
svmat ests,names(ests_tb1kwhw_)
matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_tb1kwhw_*
svmat ses,names(stderrors_tb1kwhw_)
replace stderrors_tb1kwhw_1=sqrt(stderrors_tb1kwhw_1)
capture drop CIupper_tb1kwhw_1 CIlower_tb1kwhw_1
gen CIupper_tb1kwhw_1=ests_tb1kwhw_1+stderrors_tb1kwhw_1*1.96
gen CIlower_tb1kwhw_1=ests_tb1kwhw_1-stderrors_tb1kwhw_1*1.96

tw rspike CIupper_tb1kwhw_1 CIlower_tb1kwhw_1 X_BCHfixedbins_v2 if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005,lstyle(ci) lcolor(navy) ///
	|| sc ests_tb1kwhw_1 X_BCHfixedbins_v2 if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005, mcolor(navy) lc(navy) ///
	ytitle(Mean kWh Pre-Challenge) xtitle("") xline(-0.095,lc(maroon) lp(dash)) xline(-0.1,lc(maroon) lp(solid)) ///
	xlabel(-0.2(0.05)0 0 "0%" -0.05 "-5%" -0.1 "-10%" -0.15 "-15%" -0.2 "-20%") graphregion(color(white)) ///
	legend(off) ylabel(#5) title("Pre-Challenge Use")
graph rename Balance_tb1kwhw,replace
//Figure E.3 (c)

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
eststo balance_startmonth: quietly reg start_month ibn.bins_BCHfixedbins_v2 if uniid_reduc==1 & reduc==1 & balancedset_sg6==1,nocon vce(robust)
sort X_BCHfixedbins_v2
matrix ests=e(b)'
capture drop ests_startmon_*
svmat ests,names(ests_startmon_)
matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_startmon_*
svmat ses,names(stderrors_startmon_)
replace stderrors_startmon_1=sqrt(stderrors_startmon_1)
capture drop CIupper_startmon_1 CIlower_startmon_1
gen CIupper_startmon_1=ests_startmon_1+stderrors_startmon_1*1.96
gen CIlower_startmon_1=ests_startmon_1-stderrors_startmon_1*1.96

tw rspike CIupper_startmon_1 CIlower_startmon_1 X_BCHfixedbins_v2_neg if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005,lstyle(ci) lcolor(navy) xline(0.095,lc(maroon) lp(dash)) ///
	|| con ests_startmon_1 X_BCHfixedbins_v2_neg if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005, mcolor(navy) lc(navy) ///
	ytitle(Mean Start Month) xtitle("") ///
	xlabel(0(0.05)0.2 0 "0%" 0.05 "5%" 0.1 "10%" 0.15 "15%" 0.2 "20%") graphregion(color(white)) ///
	legend(off) ylabel(#5)
graph rename BalanceStartMonth,replace
//not in paper

*-------------------------------------------------------------------------------
*Balance test of HDDs during the month or year the challenge started?
*-------------------------------------------------------------------------------

eststo balance_starthdd: quietly reg start_hdd ibn.bins_BCHfixedbins_v2 if uniid_reduc==1 & reduc==1 & balancedset_sg6==1,nocon vce(robust)
sort X_BCHfixedbins_v2
matrix ests=e(b)'
capture drop ests_starthdd_*
svmat ests,names(ests_starthdd_)
matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_starthdd_*
svmat ses,names(stderrors_starthdd_)
replace stderrors_starthdd_1=sqrt(stderrors_starthdd_1)
capture drop CIupper_starthdd_1 CIlower_starthdd_1
gen CIupper_starthdd_1=ests_starthdd_1+stderrors_starthdd_1*1.96
gen CIlower_starthdd_1=ests_starthdd_1-stderrors_starthdd_1*1.96

tw rspike CIupper_starthdd_1 CIlower_starthdd_1 X_BCHfixedbins_v2 if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005,lstyle(ci) lcolor(navy) ///
	|| sc ests_starthdd_1 X_BCHfixedbins_v2 if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005, mcolor(navy) lc(navy) ///
	ytitle(Mean HDD of Last Pre-Challenge Month) xtitle("") xline(-0.095,lc(maroon) lp(dash)) xline(-0.1,lc(maroon) lp(solid)) ///
	xlabel(-0.2(0.05)0 0 "0%" -0.05 "-5%" -0.1 "-10%" -0.15 "-15%" -0.2 "-20%") graphregion(color(white)) ///
	legend(off) ylabel(#5) title("Cold Month Before Challenge Start")
graph rename BalanceStarthdd,replace
//Figure E.3 (d)

eststo balance_startyrhdd: quietly reg start_yr_hdd ibn.bins_BCHfixedbins_v2 if uniid_reduc==1 & reduc==1 & balancedset_sg6==1,nocon vce(robust)
sort X_BCHfixedbins_v2
matrix ests=e(b)'
capture drop ests_startyrhdd_*
svmat ests,names(ests_startyrhdd_)
matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_startyrhdd_*
svmat ses,names(stderrors_startyrhdd_)
replace stderrors_startyrhdd_1=sqrt(stderrors_startyrhdd_1)
capture drop CIupper_startyrhdd_1 CIlower_startyrhdd_1
gen CIupper_startyrhdd_1=ests_startyrhdd_1+stderrors_startyrhdd_1*1.96
gen CIlower_startyrhdd_1=ests_startyrhdd_1-stderrors_startyrhdd_1*1.96

tw rspike CIupper_startyrhdd_1 CIlower_startyrhdd_1 X_BCHfixedbins_v2 if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005,lstyle(ci) lcolor(navy) ///
	|| sc ests_startyrhdd_1 X_BCHfixedbins_v2 if X_BCHfixedbins_v2>-0.205 & X_BCHfixedbins_v2<0.005, mcolor(navy) lc(navy) ///
	ytitle(Mean HDD of Pre-Program Year) xtitle("") xline(-0.095,lc(maroon) lp(dash)) xline(-0.1,lc(maroon) lp(solid)) ///
	xlabel(-0.2(0.05)0 0 "0%" -0.05 "-5%" -0.1 "-10%" -0.15 "-15%" -0.2 "-20%") graphregion(color(white)) ///
	legend(off) ylabel(#5) title("Cold Year Pre-Challenge")
graph rename BalanceStartyrhdd,replace
//not in paper
*-------------------------------------------------------------------------------
*Combine density plots
*-------------------------------------------------------------------------------
graph combine BalanceTest_Heating Balance_Floorarea Balance_tb1kwhw BalanceStarthdd,graphregion(color(white))
graph rename BalanceTests_PostCardPlot,replace
cd "$SaveOutputFiles"
graph export BalanceTests_PostCardPlot.png,replace
//Figure E.3
end

*===============================================================================
capture program drop RD_FirstStage2ndStage
program define RD_FirstStage2ndStage
*===============================================================================
*Bandwidth of 0.0025

cd "$GeneralModified"
use "temp_RD_crosssection",replace

*-------------------------------------------------------------------------------
* Generate fixed width bins of BCHchange
*-------------------------------------------------------------------------------
capture drop bins_BCHfixedbins
gen bins_BCHfixedbins=.
local binwidth=0.0025
local i=1
forvalues x=-0.095(-`binwidth')-0.30 {
	local i=`i'-1
	local x1=`x'-`binwidth'
	replace bins_BCHfixedbins=`i' if inrange(BCHchange,`x1',`x')
}
*Find the minimal value and use this to offset the set of bins
su bins_BCHfixedbins
local minval=`r(min)'
disp "`minval'"
replace bins_BCHfixedbins=bins_BCHfixedbins-`minval'+1 if bins_BCHfixedbins!=.
su bins_BCHfixedbins
local i=`r(max)'
forvalues x=-0.095(`binwidth')0.2 {
	local i=`i'+1
	local x1=`x'+`binwidth'
	replace bins_BCHfixedbins=`i' if inrange(BCHchange,`x',`x1')
}

*-------------------------------------------------------------------------------
*Save the mean values of bins_BCHfixedbins
*-------------------------------------------------------------------------------
su bins_BCHfixedbins
local maxval=`r(max)'
sort bins_BCHfixedbins

forvalues chal=1/1 {
	capture matrix drop A
	capture matrix drop temp
	tabstat BCHchange if uniid_reduc==1 & reduc==`chal' & balancedset_sg6==1,by(bins_BCHfixedbins) save
	matrix A=.
	foreach mval of numlist 1/`maxval' {
		matrix temp=r(Stat`mval')
		matrix A=A\temp
	}
	matrix A=A[2...,1]
	capture drop X_BCHfixedbins`chal'
	svmat A,names(X_BCHfixedbins`chal')
	matrix X_BCHfixedbins`chal'=A
	rename X_BCHfixedbins`chal'1 X_BCHfixedbins`chal'
}

capture drop ests_discont_BCH_*
capture drop CIupper_dc_BCH_* CIlower_dc_BCH_*
capture drop stder_dc_BCH_*

*-------------------------------------------------------------------------------
* PLOTS: of the discontinuity in continuing 
*-------------------------------------------------------------------------------

sort X_BCHfixedbins1
eststo discont_1: quietly reg Continue ibn.bins_BCHfixedbins if uniid_reduc==1 & reduc==1 & balancedset==1 & shortgapV2<=12,nocon vce(robust)
sort X_BCHfixedbins1
matrix ests=e(b)'
svmat ests,names(ests_discont_BCH_1)
rename ests_discont_BCH_11 ests_discont_BCH_1

matrix ses=e(V)'
matrix ses=vecdiag(ses)'
svmat ses,names(stder_dc_BCH_1)
rename stder_dc_BCH_11 stder_dc_BCH_1
replace stder_dc_BCH_1=sqrt(stder_dc_BCH_1)

gen CIupper_dc_BCH_1=ests_discont_BCH_1+stder_dc_BCH_1*1.96
gen CIlower_dc_BCH_1=ests_discont_BCH_1-stder_dc_BCH_1*1.96

local bs="balancedset"

eststo m13_1_both: quietly reg y2y1y0pcchange ibn.bins_BCHfixedbins if uniid_reduc==1 & reduc==1 & `bs'==1 & shortgapV2<=12,nocon vce(robust)

sort X_BCHfixedbins1
capture drop ests_m13_1_both*
matrix ests=e(b)'
svmat ests,names(ests_m13_1_both)

matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_m13_1_both*
capture drop CIupper_m13_1_both*
capture drop CIlower_m13_1_both*
svmat ses,names(stderrors_m13_1_both)
replace stderrors_m13_1_both=sqrt(stderrors_m13_1_both)
gen CIupper_m13_1_both=ests_m13_1_both+stderrors_m13_1_both*1.96
gen CIlower_m13_1_both=ests_m13_1_both-stderrors_m13_1_both*1.96

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*Plot the discontinuity in continuing
reg Continue BCHevalSuccess Center_BCHchange Center_BCHchange_s1 if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.05 & Center_BCHchange>-0.05,vce(robust)
tw rspike CIupper_dc_BCH_1 CIlower_dc_BCH_1 X_BCHfixedbins1 if X_BCHfixedbins1>-0.145 & X_BCHfixedbins1<-0.045,lstyle(ci) lcolor(navy) ///
	|| sc ests_discont_BCH_1 X_BCHfixedbins1 if X_BCHfixedbins1>-0.145 & X_BCHfixedbins1<-0.045, mcolor(navy) ///
	|| function y = _b[_cons] + _b[Center_BCHchange]*(x+0.095), range(-0.095 -0.045) lp(dash) lc(gs2) ///
	|| function y =  _b[_cons] + _b[Center_BCHchange]*(x+0.095)+ _b[BCHevalSuccess] + _b[Center_BCHchange_s1]*[x+0.095], range(-0.145 -0.095)  lp(dash) lc(gs2) ///
	ylabel(0(0.2)1,axis(1)) yscale(range(0 1)) xtitle("Credited Changes") ///
	ytitle("Probability of Re-enrolling in a 2nd Challenge",axis(1)) ///
	xlabel(-0.135(0.02)-.055 -0.135 "-13.5%" -0.115 "-11.5%" -0.095 "-9.5%" -0.075 "-7.5%" -0.055 "-5.5%") ///
	xline(-0.095,lc(maroon) lp(dash)) xline(-0.1,lc(maroon) ) graphregion(color(white)) legend(off)

graph rename Simple_S2_3_1,replace
cd "$SaveOutputFiles"
graph export Simple_S2_3_1.png,replace

local bs="balancedset"
reg y2y1y0pcchange BCHevalSuccess Center_BCHchange Center_BCHchange_s1 if reduc==1 & uniid_reduc==1 & `bs'==1 & shortgapV2<=12 & Center_BCHchange<0.05 & Center_BCHchange>-0.05,vce(robust)

tw rspike CIupper_m13_1_both CIlower_m13_1_both X_BCHfixedbins1 if X_BCHfixedbins1>-0.145 & X_BCHfixedbins1<-0.045,lstyle(ci) lcolor(navy) ///
	|| sc ests_m13_1_both X_BCHfixedbins1 if X_BCHfixedbins1>-0.145 & X_BCHfixedbins1<-0.045,mcolor(navy)  ///
	|| function y = _b[_cons] + _b[Center_BCHchange]*(x+0.095), range(-0.095 -0.045) lp(dash) lc(gs2) ///
	|| function y =  _b[_cons] + _b[Center_BCHchange]*(x+0.095)+ _b[BCHevalSuccess] + _b[Center_BCHchange_s1]*[x+0.095], range(-0.145 -0.095)  lp(dash) lc(gs2) ///
	xlabel(-0.135(0.02)-.055 -0.135 "-13.5%" -0.115 "-11.5%" -0.095 "-9.5%" -0.075 "-7.5%" -0.055 "-5.5%") ///
	xline(-0.095,lc(maroon) lp(dash)) xline(-0.1,lc(maroon)) graphregion(color(white)) legend(off) ///
	ylabel(-0.12(0.06)0.12 -0.12 "-12%" -0.06 "-6%" 0 "0%" 0.06 "6%" 0.12 "12%") ///
	ytick(-0.12(0.03)0.12) yscale(range(-0.125 0.125)) ///
	xtitle("Credited Changes") ytitle("Post-Challenge Percent Change in kWh") legend(off) graphregion(color(white)) 
	
graph rename Simple_S2_3_2,replace
cd "$SaveOutputFiles"
graph export Simple_S2_3_2.png,replace

*-------------------------------------------------------------------------------
* PLOTS: discontinuity in continue and post-challenge use alternate bandwidths
*-------------------------------------------------------------------------------
**0.0025, bw12, 0.09 window
reg Continue BCHevalSuccess Center_BCHchange Center_BCHchange_s1 if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.09 & Center_BCHchange>-0.09,vce(robust)
tw rspike CIupper_dc_BCH_1 CIlower_dc_BCH_1 X_BCHfixedbins1 if X_BCHfixedbins1>-0.185 & X_BCHfixedbins1<-0.005,lstyle(ci) lcolor(navy) ///
	|| sc ests_discont_BCH_1 X_BCHfixedbins1 if X_BCHfixedbins1>-0.185 & X_BCHfixedbins1<-0.005, mcolor(navy) ///
	|| function y = _b[_cons] + _b[Center_BCHchange]*(x+0.095), range(-0.095 -0.005) lp(dash) lc(gs2) ///
	|| function y =  _b[_cons] + _b[Center_BCHchange]*(x+0.095)+ _b[BCHevalSuccess] + _b[Center_BCHchange_s1]*[x+0.095], range(-0.185 -0.095)  lp(dash) lc(gs2) ///
	ylabel(0(0.2)1,axis(1)) yscale(range(0 1)) xtitle(Credited Changes) ///
	ytitle(Probability of Continuing to 2nd Challenge,axis(1)) ///
	xlabel(-0.175(0.04)-.015 -0.175 "-17.5%" -0.155 "-15.5%" -0.135 "-13.5%" -0.115 "-11.5%" -0.095 "-9.5%" -0.075 "-7.5%" -0.055 "-5.5%" -0.035 "-3.5%" -0.015 "-1.5%") ///
	xline(-0.095,lc(maroon) lp(dash)) xline(-0.1,lc(maroon)) graphregion(color(white)) legend(off)

graph rename Simple_S1_2_W9_bs12,replace
cd "$SaveOutputFiles"
graph export Simple_S1_2_W9_bs12.png,replace

*Window 0.09
local bs="balancedset"
reg y2y1y0pcchange BCHevalSuccess Center_BCHchange Center_BCHchange_s1 if reduc==1 & uniid_reduc==1 & `bs'==1 & shortgapV2<=12 & Center_BCHchange<0.09 & Center_BCHchange>-0.09,vce(robust)

tw rspike CIupper_m13_1_both CIlower_m13_1_both X_BCHfixedbins1 if X_BCHfixedbins1>-0.185 & X_BCHfixedbins1<-0.005 & CIupper_m13_1_both<0.3,lstyle(ci) lcolor(navy) ///
	|| sc ests_m13_1_both X_BCHfixedbins1 if X_BCHfixedbins1>-0.185 & X_BCHfixedbins1<-0.005 & CIupper_m13_1_both<0.3,mcolor(navy)  ///
	|| function y = _b[_cons] + _b[Center_BCHchange]*(x+0.095), range(-0.095 -0.005) lp(dash) lc(gs2) ///
	|| function y =  _b[_cons] + _b[Center_BCHchange]*(x+0.095)+ _b[BCHevalSuccess] + _b[Center_BCHchange_s1]*[x+0.095], range(-0.185 -0.095)  lp(dash) lc(gs2) ///
	xlabel(-0.175(0.02)-.015 -0.175 "-17.5%" -0.155 "-15.5%" -0.135 "-13.5%" -0.115 "-11.5%" -0.095 "-9.5%" -0.075 "-7.5%" -0.055 "-5.5%" -0.035 "-3.5%" -0.015 "-1.5%") ///
	xline(-0.095,lc(maroon) lp(dash)) xline(-0.1,lc(maroon)) graphregion(color(white)) legend(off) ///
	ylabel(-0.125(0.05)0.125 -0.125 "-12.5%" -0.075 "-7.5%" -0.025 "-2.5%" 0.125 "12.5%" 0.075 "7.5%" 0.025 "2.5%" ,axis(1)) ///
	xtitle("Credited Changes") ytitle("Post-Challenge Percent Change in kWh") legend(off) graphregion(color(white))

graph rename Simple_S2_2_W9_bs12,replace
cd "$SaveOutputFiles"
graph export Simple_S2_2_W9_bs12.png,replace

	
end

*===============================================================================
capture program drop RD_1stStage2ndStage_Chal23
program define RD_1stStage2ndStage_Chal23
*===============================================================================
*RF discontinuity in post-challenge use for the 2nd and 3rd challenges

cd "$GeneralModified"
use "temp_RD_crosssection",replace

*===============================================================================
* Generate fixed width bins of BCHchange
*===============================================================================
capture drop bins_BCHfixedbins
gen bins_BCHfixedbins=.
local binwidth=0.0025
local i=1
forvalues x=-0.095(-`binwidth')-0.30 {
	local i=`i'-1
	local x1=`x'-`binwidth'
	replace bins_BCHfixedbins=`i' if inrange(BCHchange,`x1',`x')
}
*Find the minimal value and use this to offset the set of bins
su bins_BCHfixedbins
local minval=`r(min)'
disp "`minval'"
replace bins_BCHfixedbins=bins_BCHfixedbins-`minval'+1 if bins_BCHfixedbins!=.
su bins_BCHfixedbins
local i=`r(max)'
forvalues x=-0.095(`binwidth')0.2 {
	local i=`i'+1
	local x1=`x'+`binwidth'
	replace bins_BCHfixedbins=`i' if inrange(BCHchange,`x',`x1')
}

*===============================================================================
* PLOTS: Post-challenge discontinuity after second challenge
*===============================================================================
*-------------------------------------------------------------------------------
*Save the mean values of bins_BCHfixedbins
*-------------------------------------------------------------------------------
su bins_BCHfixedbins
local maxval=`r(max)'
*Plot the average BCHchange by these BCHnewbins
*need to use the correct number of bins
sort bins_BCHfixedbins

forvalues chal=2/2 {
	capture matrix drop A
	capture matrix drop temp
	tabstat BCHchange if uniid_reduc==1 & reduc==`chal' & balancedset_sg6==1,by(bins_BCHfixedbins) save
	matrix A=.
	foreach mval of numlist 1/`maxval' {
		matrix temp=r(Stat`mval')
		matrix A=A\temp
	}
	matrix A=A[2...,1]
	capture drop X_BCHfixedbins`chal'
	svmat A,names(X_BCHfixedbins`chal')
	matrix X_BCHfixedbins`chal'=A
	rename X_BCHfixedbins`chal'1 X_BCHfixedbins`chal'
}

capture drop ests_discont_BCH_*
capture drop CIupper_dc_BCH_* CIlower_dc_BCH_*
capture drop stder_dc_BCH_*
*Record the P(continue) from a regression for equal width bins

*-------------------------------------------------------------------------------
*Plot for 2nd challenge
*-------------------------------------------------------------------------------

local bs="balancedset"
eststo m13_1_both: quietly reg y3y2y0pcchange ibn.bins_BCHfixedbins if uniid_reduc==1 & reduc==2 & `bs'==1 & shortgapV2<=12,nocon vce(robust)

sort X_BCHfixedbins2
capture drop ests_m13_1_both*
matrix ests=e(b)'
svmat ests,names(ests_m13_1_both)

matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_m13_1_both*
capture drop CIupper_m13_1_both*
capture drop CIlower_m13_1_both*
svmat ses,names(stderrors_m13_1_both)
replace stderrors_m13_1_both=sqrt(stderrors_m13_1_both)
gen CIupper_m13_1_both=ests_m13_1_both+stderrors_m13_1_both*1.96
gen CIlower_m13_1_both=ests_m13_1_both-stderrors_m13_1_both*1.96

local bs="balancedset"
reg y3y2y0pcchange BCHevalSuccess Center_BCHchange Center_BCHchange_s1 if reduc==2 & uniid_reduc==1 & `bs'==1 & shortgapV2<=12 & Center_BCHchange<0.05 & Center_BCHchange>-0.05,vce(robust)

tw rspike CIupper_m13_1_both CIlower_m13_1_both X_BCHfixedbins2 if X_BCHfixedbins2>-0.145 & X_BCHfixedbins2<-0.045,lstyle(ci) lcolor(navy) ///
	|| sc ests_m13_1_both X_BCHfixedbins2 if X_BCHfixedbins2>-0.145 & X_BCHfixedbins2<-0.045,mcolor(navy)  ///
	|| function y = _b[_cons] + _b[Center_BCHchange]*(x+0.095), range(-0.095 -0.045) lp(dash) lc(gs2) ///
	|| function y =  _b[_cons] + _b[Center_BCHchange]*(x+0.095)+ _b[BCHevalSuccess] + _b[Center_BCHchange_s1]*[x+0.095], range(-0.145 -0.095)  lp(dash) lc(gs2) ///
	xlabel(-0.135(0.02)-.055 -0.135 "-13.5%" -0.115 "-11.5%" -0.095 "-9.5%" -0.075 "-7.5%" -0.055 "-5.5%") ///
	xline(-0.095,lc(maroon) lp(dash)) xline(-0.1,lc(maroon)) graphregion(color(white)) legend(off) ///
	ylabel(-0.18(0.06)0.18 -0.18 "-18%" -0.12 "-12%" -0.06 "-6%" 0 "0%" 0.06 "6%" 0.12 "12%" 0.18 "18%") ///
	ytick(-0.12(0.03)0.12) yscale(range(-0.125 0.125)) ///
	xtitle("Credited Changes") ytitle("Post-Challenge Percent Change in kWh") legend(off) graphregion(color(white)) 
graph rename Simple_S3_1_1,replace
//post-challenge threshold after a second challenge

*===============================================================================
* PLOTS: Post-challenge discontinuity after third challenge
*===============================================================================
*-------------------------------------------------------------------------------
*Save the mean values of bins_BCHfixedbins
*-------------------------------------------------------------------------------
su bins_BCHfixedbins
local maxval=`r(max)'
*Plot the average BCHchange by these BCHnewbins
*need to use the correct number of bins
sort bins_BCHfixedbins

forvalues chal=3/3 {
	capture matrix drop A
	capture matrix drop temp
	tabstat BCHchange if uniid_reduc==1 & reduc==`chal' & balancedset_sg6==1,by(bins_BCHfixedbins) save
	matrix A=.
	foreach mval of numlist 1/`maxval' {
		matrix temp=r(Stat`mval')
		matrix A=A\temp
	}
	matrix A=A[2...,1]
	capture drop X_BCHfixedbins`chal'
	svmat A,names(X_BCHfixedbins`chal')
	matrix X_BCHfixedbins`chal'=A
	rename X_BCHfixedbins`chal'1 X_BCHfixedbins`chal'
}

capture drop ests_discont_BCH_*
capture drop CIupper_dc_BCH_* CIlower_dc_BCH_*
capture drop stder_dc_BCH_*

*-------------------------------------------------------------------------------
*Plot for 3rd challenge
*-------------------------------------------------------------------------------

local bs="balancedset"
eststo m13_3_both: quietly reg y4y3y0pcchange ibn.bins_BCHfixedbins if uniid_reduc==1 & reduc==3 & `bs'==1 & shortgapV2<=12,nocon vce(robust)

sort X_BCHfixedbins2
capture drop ests_m13_3_both*
matrix ests=e(b)'
svmat ests,names(ests_m13_3_both)

matrix ses=e(V)'
matrix ses=vecdiag(ses)'
capture drop stderrors_m13_3_both*
capture drop CIupper_m13_3_both*
capture drop CIlower_m13_3_both*
svmat ses,names(stderrors_m13_3_both)
replace stderrors_m13_3_both=sqrt(stderrors_m13_3_both)
gen CIupper_m13_3_both=ests_m13_3_both+stderrors_m13_3_both*1.96
gen CIlower_m13_3_both=ests_m13_3_both-stderrors_m13_3_both*1.96

local bs="balancedset"
reg y4y3y0pcchange BCHevalSuccess Center_BCHchange Center_BCHchange_s1 if reduc==3 & uniid_reduc==1 & `bs'==1 & shortgapV2<=12 & Center_BCHchange<0.05 & Center_BCHchange>-0.05,vce(robust)

tw rspike CIupper_m13_3_both CIlower_m13_3_both X_BCHfixedbins2 if X_BCHfixedbins2>-0.145 & X_BCHfixedbins2<-0.045,lstyle(ci) lcolor(navy) ///
	|| sc ests_m13_3_both X_BCHfixedbins2 if X_BCHfixedbins2>-0.145 & X_BCHfixedbins2<-0.045,mcolor(navy)  ///
	|| function y = _b[_cons] + _b[Center_BCHchange]*(x+0.095), range(-0.095 -0.045) lp(dash) lc(gs2) ///
	|| function y =  _b[_cons] + _b[Center_BCHchange]*(x+0.095)+ _b[BCHevalSuccess] + _b[Center_BCHchange_s1]*[x+0.095], range(-0.145 -0.095)  lp(dash) lc(gs2) ///
	xlabel(-0.135(0.02)-.055 -0.135 "-13.5%" -0.115 "-11.5%" -0.095 "-9.5%" -0.075 "-7.5%" -0.055 "-5.5%") ///
	xline(-0.095,lc(maroon) lp(dash)) xline(-0.1,lc(maroon)) graphregion(color(white)) legend(off) ///
	ylabel(-0.18(0.06)0.24 -0.18 "-18%" -0.12 "-12%" -0.06 "-6%" 0 "0%" 0.06 "6%" 0.12 "12%" 0.18 "18%" 0.24 "24%") ///
	ytick(-0.12(0.03)0.12) yscale(range(-0.125 0.125)) ///
	xtitle("Credited Changes") ytitle("Post-Challenge Percent Change in kWh") legend(off) graphregion(color(white)) 
graph rename Simple_S3_3_1,replace
//post-challenge threshold after a second challenge

end

*===============================================================================
capture program drop RD_McCrary
program define RD_McCrary
*===============================================================================

cd "$GeneralModified"
use "temp_RD_crosssection",replace

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
* Histograms
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*Histogram using balancedset_change==1
hist BCHchange if uniid_reduc==1 & reduc==1 & heating!=3 & BCHchange>=-0.2 & BCHchange<=-0.0 & balancedset_change==1, ///
	xlabel(-0.2(0.05)0 -0.2 "-20%" -0.15 "-15%" -0.1 "-10%" -0.05 "-5%" 0 "0%") xtitle("") /// 
	xline(-0.095,lc(maroon) lp(dash)) xline(-0.1,lc(maroon)) w(0.0025)  graphregion(color(white)) fcolor(gs10) lcolor(gs14)

graph rename McCrary_B_5_4_3,replace

*-------------------------------------------------------------------------------
*Households target 10%: Use that -0.1 as the cutoff
*-------------------------------------------------------------------------------

capture drop Xj Yj r0 fhat se_fhat
DCdensity BCHchange if BCHchange>-0.2 & BCHchange<0.0 & balancedset_change==1, breakpoint(-0.095) generate(Xj Yj r0 fhat se_fhat)
												 
capture drop fhat_*
gen fhat_upper=fhat+1.96*se_fhat
gen fhat_lower=fhat-1.96*se_fhat

tw line fhat r0 if r0<-0.095 & r0>-0.2,lc(black) ///
	|| line fhat r0 if r0>-0.095 & r0<0,lc(black) ///
	|| rline fhat_upper fhat_lower r0 if r0<-0.095 & r0>-0.2,lc(gs4) lw(vthin) ///
	|| rline fhat_upper fhat_lower r0 if r0>-0.095 & r0<0,lc(gs4) lw(vthin) ///
	|| sc Yj Xj if Xj>-0.2 & Xj<0,m(circle_hollow) mc(gs10) msize(small)  ///
	xlabel(-0.2(0.05)0 -0.2 "-20%" -0.15 "-15%" -0.1 "-10%" -0.05 "-5%" 0 "0%") ///
	xtitle("Percent Change in Annual kWh During First Challenge (Credited Changes)") ytitle("Density") ///
	legend(off) graphregion(color(white)) xline(-0.095,lc(maroon) lp(dash) lw(thin)) xline(-0.1) 

graph rename McCrary_BCHchange_B_4_5,replace
cd "$SaveOutputFiles"
graph export McCrary_BCHchange_B_4_5.png,replace

*For 10% threshold
capture drop Xj Yj r0 fhat se_fhat
DCdensity BCHchange if BCHchange>-0.2 & BCHchange<0.0 & balancedset_change==1, breakpoint(-0.1) generate(Xj Yj r0 fhat se_fhat)

							 
capture drop fhat_*
gen fhat_upper=fhat+1.96*se_fhat
gen fhat_lower=fhat-1.96*se_fhat

tw line fhat r0 if r0<-0.1 & r0>-0.2,lc(black) ///
	|| line fhat r0 if r0>-0.1 & r0<0,lc(black) ///
	|| rline fhat_upper fhat_lower r0 if r0<-0.1 & r0>-0.2,lc(gs4) lw(vthin) ///
	|| rline fhat_upper fhat_lower r0 if r0>-0.1 & r0<0,lc(gs4) lw(vthin) ///
	|| sc Yj Xj if Xj>-0.2 & Xj<0,m(circle_hollow) mc(gs10) msize(small)  ///
	xlabel(-0.2(0.05)0 -0.2 "-20%" -0.15 "-15%" -0.1 "-10%" -0.05 "-5%" 0 "0%") ///
	xtitle("Percent Change in Annual kWh During First Challenge (Credited Changes)") ytitle("Density") ///
	legend(off) graphregion(color(white)) xline(-0.095,lc(maroon) lp(dash) lw(thin)) xline(-0.1) 

graph rename McCrary_BCHchange_B_4_6,replace

graph combine McCrary_B_5_4_3 McCrary_BCHchange_B_4_6,r(2) graphregion(color(white)) iscale(0.75) imargin(vsmall) name(McCrary_DensityCheck,replace)
graph display McCrary_DensityCheck, xsize(9) ysize(10)
cd "$SaveOutputFiles"
graph export McCrary_DensityCheck.png,replace

end

*===============================================================================
capture program drop RD_Ests
program define RD_Ests
*===============================================================================

cd "$GeneralModified"
use "temp_RD_crosssection",replace

*-------------------------------------------------------------------------------
* Table 3
*-------------------------------------------------------------------------------
eststo ivreg_y2y1y0_m3_3: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.03 & Center_BCHchange>-0.03,savefprefix(ivy2y1y0_m33) robust
eststo ivreg_y2y1y0_m3_4: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.04 & Center_BCHchange>-0.04,savefprefix(ivy2y1y0_m34) robust
eststo ivreg_y2y1y0_m3_5: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.05 & Center_BCHchange>-0.05,savefprefix(ivy2y1y0_m35) robust saverf
eststo ivreg_y2y1y0_m3_6: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.06 & Center_BCHchange>-0.06,savefprefix(ivy2y1y0_m36) robust
eststo ivreg_y2y1y0_m3_7: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.07 & Center_BCHchange>-0.07,savefprefix(ivy2y1y0_m37) robust

estadd local bandwidth "3pc" , replace : ivreg_y2y1y0_m3_3
estadd local bandwidth "4pc" , replace : ivreg_y2y1y0_m3_4
estadd local bandwidth "5pc" , replace : ivreg_y2y1y0_m3_5
estadd local bandwidth "6pc" , replace : ivreg_y2y1y0_m3_6
estadd local bandwidth "7pc" , replace : ivreg_y2y1y0_m3_7

foreach regname in 3 4 5 6 7 {
	disp("`regname'")
	estimates restore ivreg_y2y1y0_m3_`regname'
	*mat list e(first)
	mat first=e(first)
	local fsstat=first[4,1]
	disp "`fsstat'"
	estadd scalar fs=`fsstat', replace: ivreg_y2y1y0_m3_`regname' ivy2y1y0_m3`regname'*
}


*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*Second stage
esttab ivreg_y2y1y0_m3_3 ivreg_y2y1y0_m3_4 ivreg_y2y1y0_m3_5 ivreg_y2y1y0_m3_6 ivreg_y2y1y0_m3_7, compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(bandwidth fs N, label("Window" "F-stat"))
//Table 3 Panel B

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*First stage

*For tex output, need some blank columns. Also want to start with 9% window and go smaller
gen y2y1y0fake = y2y1y0pcchange*0.001
eststo Blank1: quietly reg y2y1y0pcchange y2y1y0fake if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.05 & Center_BCHchange>-0.05,nocon vce(robust)

esttab Blank1 Blank1 ivy2y1y0_m37Continue ivy2y1y0_m36Continue ivy2y1y0_m35Continue ivy2y1y0_m34Continue ivy2y1y0_m33Continue,drop(y2y1y0fake) compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(fs, label( "F-stat")) 
//Table 3 Panel A

esttab Blank1 Blank1 ivy2y1y0_m37Continue ivy2y1y0_m36Continue ivy2y1y0_m35Continue ivy2y1y0_m34Continue ivy2y1y0_m33Continue,drop(y2y1y0fake) compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(fs, label( "F-stat"))  tex fragment
//Table 3 Panel A

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*Include OLS and FS results that correspond to this specification
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
eststo OLSm3_BCH_a: quietly reg y2y1y0pcchange Continue if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12, vce(robust)
eststo OLSm3_BCH_w5: quietly reg y2y1y0pcchange Continue Center_BCHchange Center_BCHchange_s1 Center_mychange if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.05 & Center_BCHchange>-0.05, vce(robust)

estadd local bandwidth "5pc" , replace : OLSm3_BCH_w5

*include RF instead of 5% OLS
esttab OLSm3_BCH_a _ivreg2_y2y1y0pcchange ivreg_y2y1y0_m3_7 ivreg_y2y1y0_m3_6 ivreg_y2y1y0_m3_5 ivreg_y2y1y0_m3_4 ivreg_y2y1y0_m3_3, compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(bandwidth fs N, label("Window" "F-stat"))
//Table 3 Panel B

esttab OLSm3_BCH_a _ivreg2_y2y1y0pcchange ivreg_y2y1y0_m3_7 ivreg_y2y1y0_m3_6 ivreg_y2y1y0_m3_5 ivreg_y2y1y0_m3_4 ivreg_y2y1y0_m3_3, compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(bandwidth fs N, label("Window" "F-stat")) tex fragment
//Table 3 Panel B

*-------------------------------------------------------------------------------
*Robustness checks: include other covariates
*-------------------------------------------------------------------------------
capture frame drop hdd
frame create hdd
frame change hdd
cd "$datafrom"
use hdd_data,clear

*Generate rolling annual changes in hdd defined by the start of a challenge, then merge these in with the rd observations
capture drop cum*
capture drop id_tot*
capture drop temp*

tsset year_month
sort year_month
gen cum_sum=sum(hdd)

gen id_tot_hdd_b1=L1.cum_sum-L13.cum_sum
replace id_tot_hdd_b1=L1.cum_sum if year_month==ym(2006,1)
gen id_tot_hdd_a0=F11.cum_sum-L1.cum_sum
gen id_tot_hdd_a1=F23.cum_sum-F11.cum_sum

gen hddchange_tb1ta0=(id_tot_hdd_a0-id_tot_hdd_b1)/id_tot_hdd_b1
gen hddchange_ta0ta1=(id_tot_hdd_a1-id_tot_hdd_a0)/id_tot_hdd_a0

frame change default
*want the percent change in hdd for the 12 months before starting the second challenge. The gap length makes this hard to calculate.
rename begin_challenge_1 year_month
frlink m:1 year_month,frame(hdd)
frget hddchange_tb1ta0=hddchange_tb1ta0,from(hdd)
rename year_month begin_challenge_1
drop hdd

*Verify that I can produce the main estimates correctly
eststo test_y2y1y0_m3_7: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.07 & Center_BCHchange>-0.07,savefprefix(ivy2y1y0_m37) robust
*2nd stage coef on Continue is -.126 for 7% bandwidth. This is very close to the main estimates

*Check 5% BW on main estimates
eststo test_y2y1y0_m3_5: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.05 & Center_BCHchange>-0.05,savefprefix(ivy2y1y0_m35) robust saverf
*Coef on continue is -0.23, which matches the paper estimates

*include covariates in the estimation
eststo ivreg_y2y1y0_m4_3: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange i.heating i.buildingtype hddchange_tb1ta0 (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.03 & Center_BCHchange>-0.03,savefprefix(ivy2y1y0_m43) robust
eststo ivreg_y2y1y0_m4_4: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange i.heating i.buildingtype hddchange_tb1ta0 (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.04 & Center_BCHchange>-0.04,savefprefix(ivy2y1y0_m44) robust
eststo ivreg_y2y1y0_m4_5: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange i.heating i.buildingtype hddchange_tb1ta0 (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.05 & Center_BCHchange>-0.05,savefprefix(ivy2y1y0_m45) robust saverf
eststo ivreg_y2y1y0_m4_6: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange i.heating i.buildingtype hddchange_tb1ta0 (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.06 & Center_BCHchange>-0.06,savefprefix(ivy2y1y0_m46) robust
eststo ivreg_y2y1y0_m4_7: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange i.heating i.buildingtype hddchange_tb1ta0 (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.07 & Center_BCHchange>-0.07,savefprefix(ivy2y1y0_m47) robust

estadd local bandwidth "3pc" , replace : ivreg_y2y1y0_m4_3
estadd local bandwidth "4pc" , replace : ivreg_y2y1y0_m4_4
estadd local bandwidth "5pc" , replace : ivreg_y2y1y0_m4_5
estadd local bandwidth "6pc" , replace : ivreg_y2y1y0_m4_6
estadd local bandwidth "7pc" , replace : ivreg_y2y1y0_m4_7

foreach regname in 3 4 5 6 7 {
	disp("`regname'")
	estimates restore ivreg_y2y1y0_m4_`regname'
	*mat list e(first)
	mat first=e(first)
	local fsstat=first[4,1]
	disp "`fsstat'"
	estadd scalar fs=`fsstat', replace: ivreg_y2y1y0_m4_`regname' ivy2y1y0_m4`regname'*
}


*Additional covariates
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*Second stage
esttab ivreg_y2y1y0_m4_7 ivreg_y2y1y0_m4_6 ivreg_y2y1y0_m4_5 ivreg_y2y1y0_m4_4 ivreg_y2y1y0_m4_3, compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(bandwidth fs N, label("Window" "F-stat")) keep(Continue Center* hdd* _cons)

esttab ivreg_y2y1y0_m4_7 ivreg_y2y1y0_m4_6 ivreg_y2y1y0_m4_5 ivreg_y2y1y0_m4_4 ivreg_y2y1y0_m4_3, compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(bandwidth fs N, label("Window" "F-stat")) keep(Continue Center* hdd* _cons) tex fragment
//Table 6.1, Panel B

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*First stage
esttab ivy2y1y0_m47Continue ivy2y1y0_m46Continue ivy2y1y0_m45Continue ivy2y1y0_m44Continue ivy2y1y0_m43Continue,compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(fs, label( "F-stat")) keep(BCHeval* Center* hdd* _cons)
//Table 6.1, Panel B

*tex output
esttab ivy2y1y0_m47Continue ivy2y1y0_m46Continue ivy2y1y0_m45Continue ivy2y1y0_m44Continue ivy2y1y0_m43Continue,compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(fs, label( "F-stat")) keep(BCHeval* Center* hdd* _cons) tex
//Table 6.1, Panel B

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*Restricted billing
eststo ivreg_y2y1y0_m5_3: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.03 & Center_BCHchange>-0.03 & Center_mychange<0.05 & Center_mychange>-0.05 ,savefprefix(ivy2y1y0_m53) robust
eststo ivreg_y2y1y0_m5_4: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.04 & Center_BCHchange>-0.04 & Center_mychange<0.05 & Center_mychange>-0.05,savefprefix(ivy2y1y0_m54) robust
eststo ivreg_y2y1y0_m5_5: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.05 & Center_BCHchange>-0.05 & Center_mychange<0.05 & Center_mychange>-0.05,savefprefix(ivy2y1y0_m55) robust saverf
eststo ivreg_y2y1y0_m5_6: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.06 & Center_BCHchange>-0.06 & Center_mychange<0.05 & Center_mychange>-0.05,savefprefix(ivy2y1y0_m56) robust
eststo ivreg_y2y1y0_m5_7: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=12 & Center_BCHchange<0.07 & Center_BCHchange>-0.07 & Center_mychange<0.05 & Center_mychange>-0.05,savefprefix(ivy2y1y0_m57) robust

estadd local bandwidth "3pc" , replace : ivreg_y2y1y0_m5_3
estadd local bandwidth "4pc" , replace : ivreg_y2y1y0_m5_4
estadd local bandwidth "5pc" , replace : ivreg_y2y1y0_m5_5
estadd local bandwidth "6pc" , replace : ivreg_y2y1y0_m5_6
estadd local bandwidth "7pc" , replace : ivreg_y2y1y0_m5_7

foreach regname in 3 4 5 6 7 {
	disp("`regname'")
	estimates restore ivreg_y2y1y0_m5_`regname'
	*mat list e(first)
	mat first=e(first)
	local fsstat=first[4,1]
	disp "`fsstat'"
	estadd scalar fs=`fsstat', replace: ivreg_y2y1y0_m5_`regname' ivy2y1y0_m5`regname'*
}


*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*Second stage
esttab ivreg_y2y1y0_m5_7 ivreg_y2y1y0_m5_6 ivreg_y2y1y0_m5_5 ivreg_y2y1y0_m5_4 ivreg_y2y1y0_m5_3, compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(bandwidth fs N, label("Window" "F-stat")) 

esttab ivreg_y2y1y0_m5_7 ivreg_y2y1y0_m5_6 ivreg_y2y1y0_m5_5 ivreg_y2y1y0_m5_4 ivreg_y2y1y0_m5_3, compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(bandwidth fs N, label("Window" "F-stat")) tex
//Table G.2 Panel B

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*First stage
esttab ivy2y1y0_m57Continue ivy2y1y0_m56Continue ivy2y1y0_m55Continue ivy2y1y0_m54Continue ivy2y1y0_m53Continue,compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(fs, label( "F-stat"))

esttab ivy2y1y0_m57Continue ivy2y1y0_m56Continue ivy2y1y0_m55Continue ivy2y1y0_m54Continue ivy2y1y0_m53Continue,compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(fs, label( "F-stat")) tex
//Table G.2 Panel A


*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*1st order Bias-corrected estimates
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) h(0.03) b(0.03) all
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) h(0.04) b(0.04) all
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) h(0.05) b(0.05) all
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) h(0.06) b(0.06) all
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) h(0.07) b(0.07) all
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) h(0.08) b(0.08) all
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) h(0.09) b(0.09) all
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) all

*2nd order bias corrected
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) h(0.03) b(0.03) all p(2) q(3)
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) h(0.04) b(0.04) all p(2) q(3)
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) h(0.05) b(0.05) all p(2) q(3)
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) h(0.06) b(0.06) all p(2) q(3)
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) h(0.07) b(0.07) all p(2) q(3)
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) h(0.08) b(0.08) all p(2) q(3)
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) h(0.09) b(0.09) all p(2) q(3)
rdrobust y2y1y0pcchange Center_BCHchange if balancedset==1 & shortgapV2<=12,c(0) fuzzy(Continue) all

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*Estimates using the 6 month gap

eststo ivreg_y2y1y0_m6_3: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=6 & Center_BCHchange<0.03 & Center_BCHchange>-0.03,savefprefix(ivy2y1y0_m63) robust
eststo ivreg_y2y1y0_m6_4: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=6 & Center_BCHchange<0.04 & Center_BCHchange>-0.04,savefprefix(ivy2y1y0_m64) robust
eststo ivreg_y2y1y0_m6_5: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=6 & Center_BCHchange<0.05 & Center_BCHchange>-0.05,savefprefix(ivy2y1y0_m65) robust saverf
eststo ivreg_y2y1y0_m6_6: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=6 & Center_BCHchange<0.06 & Center_BCHchange>-0.06,savefprefix(ivy2y1y0_m66) robust
eststo ivreg_y2y1y0_m6_7: ivreg2 y2y1y0pcchange Center_BCHchange Center_BCHchange_s1 Center_mychange (Continue=BCHevalSuccess) if reduc==1 & uniid_reduc==1 & balancedset==1 & shortgapV2<=6 & Center_BCHchange<0.07 & Center_BCHchange>-0.07,savefprefix(ivy2y1y0_m67) robust

estadd local bandwidth "3pc" , replace : ivreg_y2y1y0_m6_3
estadd local bandwidth "4pc" , replace : ivreg_y2y1y0_m6_4
estadd local bandwidth "5pc" , replace : ivreg_y2y1y0_m6_5
estadd local bandwidth "6pc" , replace : ivreg_y2y1y0_m6_6
estadd local bandwidth "7pc" , replace : ivreg_y2y1y0_m6_7

foreach regname in 3 4 5 6 7 {
	disp("`regname'")
	estimates restore ivreg_y2y1y0_m6_`regname'
	*mat list e(first)
	mat first=e(first)
	local fsstat=first[4,1]
	disp "`fsstat'"
	estadd scalar fs=`fsstat', replace: ivreg_y2y1y0_m6_`regname' ivy2y1y0_m6`regname'*
}


*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*Second stage
esttab ivreg_y2y1y0_m6_7 ivreg_y2y1y0_m6_6 ivreg_y2y1y0_m6_5 ivreg_y2y1y0_m6_4 ivreg_y2y1y0_m6_3, compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(bandwidth fs N, label("Window" "F-stat")) tex
//Table G.5 - Panel B

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*First stage
esttab ivy2y1y0_m67Continue ivy2y1y0_m66Continue ivy2y1y0_m65Continue ivy2y1y0_m64Continue ivy2y1y0_m63Continue,compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(fs, label( "F-stat")) tex
//Table G.5 - Panel A

end

*===============================================================================
capture program drop iv_lnkwhw
program define iv_lnkwhw
*===============================================================================
* Estimating using lnkwhw specification

cd "$GeneralModified"
use "BalancedData",replace

capture drop _est*
*-------------------------------------------------------------------------------
* Simplify dataset from the full generate data program
*-------------------------------------------------------------------------------
*drop indicators and merge in the total_kwh_num and similar values
capture drop t_t*
foreach x of numlist 0(1)10 {
	capture drop t_b`x'_fm*
	capture drop t_a`x'_fm*
}
foreach x of numlist 3(1)10 {
	capture drop t_b`x'*
	capture drop t_a`x'*
}

drop reduc_gap_*
drop lastmonth
drop reward

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
* Generate mychange
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
gen prepostyears=0 if t_a0==1
replace prepostyears=-1 if t_b1==1
replace prepostyears=1 if t_a1==1
replace prepostyears=2 if t_a2==1

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
* Generate prepostyears_v2
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*need the 12 months before a second challenge, or 12 months after the 1st challenge
tsset id date
gen prepostyears_v2=1 if t_a1==1
*Assign prepostyears_v2=0 if 12 months before a second challenge
replace prepostyears_v2=0 if F12.t_a1==1
*Assign prepostyears_v2=0 if 12 months after a 1st challenge and there is no second challenge
replace prepostyears_v2=0 if L12.t_a0==1 & num_challenges==1
*br id tps reduc date t_b2 t_b1 t_a0 t_a1 prepost* if id==20003


*bysort id prepostyears: gen uniprepost=1 if _n==1
bysort id prepostyears: gen num_months=_N if prepostyears!=.

capture drop temp
forvalues i=0/1 {
	disp "`i'"
	capture drop temp
	bysort id: egen temp=total(kwhw) if t_a`i'==1 & num_months==12
	bysort id: egen ta`i'_year_kwhw_1=mean(temp)
}
forvalues i=1/1 {
	disp "`i'"
	capture drop temp
	bysort id: egen temp=total(kwhw) if t_b`i'==1 & num_months==12
	bysort id: egen tb`i'_year_kwhw_1=mean(temp)
}


gen mychange=(ta0_year_kwhw_1-tb1_year_kwhw_1)/(tb1_year_kwhw_1) if reduc==1
bysort id: egen id_mychange=mean(mychange)

gen Center_mychange=mychange+0.095

*-------------------------------------------------------------------------------
* Indicator for whether HH continues past the current id or not
*-------------------------------------------------------------------------------
gen Continue=1 if num_challenges>reduc & reduc!=0
replace Continue=0 if Continue==.

bysort id reduc: gen uniid_reduc=1 if _n==1

*Merge in info needed like total_kwh_num
rename reduc challenge_order
merge m:1 id challenge_order using Accounts_All,keepusing(total_kwh_num adjusted_historical_kwh_num adjusted_historical_kwh_num success2) keep(1 3) nogen
rename challenge_order reduc

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*Windows around threshold
gen BCHchange=(total_kwh_num-adjusted_historical_kwh_num)/adjusted_historical_kwh_num if reduc!=0
gen BCHevalSuccess=1 if BCHchange<-0.095
replace BCHevalSuccess=0 if BCHevalSuccess==.

tab BCHevalSuccess success2

capture drop temp
capture drop id_BCHchange_1
bysort id: gen temp=BCHchange if reduc==1
bysort id: egen id_BCHchange_1=mean(temp)
drop temp
su BCHchange if id_BCHchange_1>-0.15 & id_BCHchange_1<-0.05 & reduc==1


gen Center_id_BCHchange_1=id_BCHchange_1+0.095

*Generate indicators needed
gen Continue_ta0_ta1=t_a1 if num_challenges>1
replace Continue_ta0_ta1=0 if Continue_ta0_ta1==.

capture drop temp*
capture drop id_BCHSuccess_t_a1
bysort id: gen temp=BCHevalSuccess if t_a0==1
bysort id: egen temp2=mean(temp)
bysort id: gen id_BCHSuccess_t_a1=temp2 if t_a1==1
drop temp
replace id_BCHSuccess_t_a1=0 if id_BCHSuccess_t_a1==.

gen id_BCHchange_t_a1=id_BCHchange_1 if t_a1==1
replace id_BCHchange_t_a1=0 if id_BCHchange_t_a1==.

gen id_BCHchange_cut=id_BCHchange_1+0.095

gen id_BCHchange_cut_t_a1=id_BCHchange_cut if t_a1==1
replace id_BCHchange_cut_t_a1=0 if id_BCHchange_cut_t_a1==.

gen id_BCHchange_cut_s1=id_BCHchange_cut if BCHevalSuccess==1
replace id_BCHchange_cut_s1=0 if id_BCHchange_cut_s1==.

gen id_BCHchange_cut_t_a1_s1=id_BCHchange_cut_t_a1 if BCHevalSuccess==1
replace id_BCHchange_cut_t_a1_s1=0 if  id_BCHchange_cut_t_a1_s1==.

*===============================================================================
* ivreg using lnkwhw
*===============================================================================

eststo iv_lnkwhw_3: xtivreg2 lnkwhw t_a1 id_BCHchange_cut_t_a1 id_BCHchange_cut_t_a1_s1 (Continue_ta0_ta1 = id_BCHSuccess_t_a1) if inlist(prepostyears_v2,0,1) & Center_id_BCHchange_1>-0.03 & Center_id_BCHchange_1<0.03 & balancedset==1 & shortgapV2<=12,fe cl(id) savefprefix(iv_ln_fs_3)

eststo iv_lnkwhw_4: xtivreg2 lnkwhw t_a1 id_BCHchange_cut_t_a1 id_BCHchange_cut_t_a1_s1 (Continue_ta0_ta1 = id_BCHSuccess_t_a1) if inlist(prepostyears_v2,0,1) & Center_id_BCHchange_1>-0.04 & Center_id_BCHchange_1<0.04 & balancedset==1 & shortgapV2<=12,fe cl(id) savefprefix(iv_ln_fs_4)

eststo iv_lnkwhw_5: xtivreg2 lnkwhw t_a1 id_BCHchange_cut_t_a1 id_BCHchange_cut_t_a1_s1 (Continue_ta0_ta1 = id_BCHSuccess_t_a1) if inlist(prepostyears_v2,0,1) & Center_id_BCHchange_1>-0.05 & Center_id_BCHchange_1<0.05 & balancedset==1 & shortgapV2<=12,fe cl(id) savefprefix(iv_ln_fs_5)

eststo iv_lnkwhw_6: xtivreg2 lnkwhw t_a1 id_BCHchange_cut_t_a1 id_BCHchange_cut_t_a1_s1 (Continue_ta0_ta1 = id_BCHSuccess_t_a1) if inlist(prepostyears_v2,0,1) & Center_id_BCHchange_1>-0.06 & Center_id_BCHchange_1<0.06 & balancedset==1 & shortgapV2<=12,fe cl(id) savefprefix(iv_ln_fs_6)

eststo iv_lnkwhw_7: xtivreg2 lnkwhw t_a1 id_BCHchange_cut_t_a1 id_BCHchange_cut_t_a1_s1 (Continue_ta0_ta1 = id_BCHSuccess_t_a1) if inlist(prepostyears_v2,0,1) & Center_id_BCHchange_1>-0.07 & Center_id_BCHchange_1<0.07 & balancedset==1 & shortgapV2<=12,fe cl(id) savefprefix(iv_ln_fs_7)

foreach regname in 7 6 5 4 3 {
	disp("`regname'")
	estimates restore iv_lnkwhw_`regname'
	*mat list e(first)
	mat first=e(first)
	local fsstat=first[4,1]
	disp "`fsstat'"
	estadd scalar fs=`fsstat', replace: iv_lnkwhw_`regname' iv_ln_fs_`regname'*
	local test= e(N_clust) 
	estadd scalar nids=`test', replace: iv_lnkwhw_`regname' iv_ln_fs_`regname'*
}

esttab iv_lnkwhw_7 iv_lnkwhw_6 iv_lnkwhw_5 iv_lnkwhw_4 iv_lnkwhw_3, order(Continue_ta0_ta1 t_a1) compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(N nids)

esttab iv_ln_fs_7Continue_ta0_ta1 iv_ln_fs_6Continue_ta0_ta1 iv_ln_fs_5Continue_ta0_ta1 iv_ln_fs_4Continue_ta0_ta1 iv_ln_fs_3Continue_ta0_ta1, compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(fs, label( "F-stat")) 

*OLS for plot
eststo OLS_lnk_m3_a: xtivreg2 lnkwhw t_a1 Continue_ta0_ta1 if inlist(prepostyears_v2,0,1) & balancedset==1 & shortgapV2<=12,fe cl(id)
eststo OLS_lnk_m3_w5: xtivreg2 lnkwhw t_a1 Continue_ta0_ta1 if inlist(prepostyears_v2,0,1) & Center_id_BCHchange_1>-0.05 & Center_id_BCHchange_1<0.05 & balancedset==1 & shortgapV2<=12,fe cl(id)

*Blank indicators
eststo lnk_Blanks1: xtivreg2 lnkwhw t_a2 if inlist(prepostyears,0,1,2) & balancedset==1 & shortgapV2<=12,fe cl(id)

*2nd stage
esttab OLS_lnk_m3_a OLS_lnk_m3_w5 lnk_Blanks1 iv_lnkwhw_7 iv_lnkwhw_6 iv_lnkwhw_5 iv_lnkwhw_4 iv_lnkwhw_3, order(Continue_ta0_ta1 t_a1) compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(N nids) tex
//Table G.6 Panel B

*1st stage
esttab iv_ln_fs_7Continue_ta0_ta1 iv_ln_fs_6Continue_ta0_ta1 iv_ln_fs_5Continue_ta0_ta1 iv_ln_fs_4Continue_ta0_ta1 iv_ln_fs_3Continue_ta0_ta1, compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(fs, label( "F-stat")) tex
//Table G.6 Panel A

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*6 month gap
eststo iv_lnkwhw6_3: xtivreg2 lnkwhw t_a1 id_BCHchange_cut_t_a1 id_BCHchange_cut_t_a1_s1 (Continue_ta0_ta1 = id_BCHSuccess_t_a1) if inlist(prepostyears_v2,0,1) & Center_id_BCHchange_1>-0.03 & Center_id_BCHchange_1<0.03 & balancedset==1 & shortgapV2<=6,fe cl(id) savefprefix(iv_fs6_3)

eststo iv_lnkwhw6_4: xtivreg2 lnkwhw t_a1 id_BCHchange_cut_t_a1 id_BCHchange_cut_t_a1_s1 (Continue_ta0_ta1 = id_BCHSuccess_t_a1) if inlist(prepostyears_v2,0,1) & Center_id_BCHchange_1>-0.04 & Center_id_BCHchange_1<0.04 & balancedset==1 & shortgapV2<=6,fe cl(id) savefprefix(iv_fs6_4)

eststo iv_lnkwhw6_5: xtivreg2 lnkwhw t_a1 id_BCHchange_cut_t_a1 id_BCHchange_cut_t_a1_s1 (Continue_ta0_ta1 = id_BCHSuccess_t_a1) if inlist(prepostyears_v2,0,1) & Center_id_BCHchange_1>-0.05 & Center_id_BCHchange_1<0.05 & balancedset==1 & shortgapV2<=6,fe cl(id) savefprefix(iv_fs6_5)

eststo iv_lnkwhw6_6: xtivreg2 lnkwhw t_a1 id_BCHchange_cut_t_a1 id_BCHchange_cut_t_a1_s1 (Continue_ta0_ta1 = id_BCHSuccess_t_a1) if inlist(prepostyears_v2,0,1) & Center_id_BCHchange_1>-0.06 & Center_id_BCHchange_1<0.06 & balancedset==1 & shortgapV2<=6,fe cl(id) savefprefix(iv_fs6_6)

eststo iv_lnkwhw6_7: xtivreg2 lnkwhw t_a1 id_BCHchange_cut_t_a1 id_BCHchange_cut_t_a1_s1 (Continue_ta0_ta1 = id_BCHSuccess_t_a1) if inlist(prepostyears_v2,0,1) & Center_id_BCHchange_1>-0.07 & Center_id_BCHchange_1<0.07 & balancedset==1 & shortgapV2<=6,fe cl(id) savefprefix(iv_fs6_7)

foreach regname in 7 6 5 4 3 {
	disp("`regname'")
	estimates restore iv_lnkwhw6_`regname'
	*mat list e(first)
	mat first=e(first)
	local fsstat=first[4,1]
	disp "`fsstat'"
	estadd scalar fs=`fsstat', replace: iv_lnkwhw6_`regname' iv_fs6_`regname'*
	local test= e(N_clust) 
	estadd scalar nids=`test', replace: iv_lnkwhw6_`regname' iv_fs6_`regname'*
}

esttab iv_lnkwhw6_7 iv_lnkwhw6_6 iv_lnkwhw6_5 iv_lnkwhw6_4 iv_lnkwhw6_3, order(Continue_ta0_ta1 t_a1) compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(N nids)

esttab iv_fs6_7Continue_ta0_ta1 iv_fs6_6Continue_ta0_ta1 iv_fs6_5Continue_ta0_ta1 iv_fs6_4Continue_ta0_ta1 iv_fs6_3Continue_ta0_ta1, compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(fs, label( "F-stat")) 

*OLS for plot
eststo OLS_lnk_m3_a: xtivreg2 lnkwhw t_a1 Continue_ta0_ta1 if inlist(prepostyears_v2,0,1) & balancedset==1 & shortgapV2<=6,fe cl(id)
eststo OLS_lnk_m3_w5: xtivreg2 lnkwhw t_a1 Continue_ta0_ta1 if inlist(prepostyears_v2,0,1) & Center_id_BCHchange_1>-0.05 & Center_id_BCHchange_1<0.05 & balancedset==1 & shortgapV2<=6,fe cl(id)

*Blank indicators
eststo lnk_Blanks1: xtivreg2 lnkwhw t_a2 if inlist(prepostyears,0,1,2) & balancedset==1 & shortgapV2<=6,fe cl(id)

*2nd stage
esttab OLS_lnk_m3_a OLS_lnk_m3_w5 lnk_Blanks1 iv_lnkwhw6_7 iv_lnkwhw6_6 iv_lnkwhw6_5 iv_lnkwhw6_4 iv_lnkwhw6_3, order(Continue_ta0_ta1 t_a1) compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(N nids) tex
//Table G.7 Panel B

*1st stage
esttab iv_fs6_7Continue_ta0_ta1 iv_fs6_6Continue_ta0_ta1 iv_fs6_5Continue_ta0_ta1 iv_fs6_4Continue_ta0_ta1 iv_fs6_3Continue_ta0_ta1, compress nogaps star(* 0.10 ** 0.05 *** 0.01) se mti s(fs, label( "F-stat")) tex
//Table G.7 Panel A

end

*========================================================================
*Call the main program to run this file.
*========================================================================
main
